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ABSTRACT 


This  paper  presents  a  survey  of  the  current  state-of-the-art  in  multidisciplinary  aeromechanical  analyses 
which  integrate  advanced  Computational  Structural  Dynamics  (CSD)  and  Computational  Fluid  Dynamics 
(CFD)  methods.  The  application  areas  to  be  surveyed  include  fixed  wing  aircraft,  turbomachinary,  and 
rotary  wing  aircraft.  The  objective  of  the  authors  in  the  present  paper  —  together  with  a  companion  paper 
on  requirements  —  is  to  lay  out  a  path  for  a  High  Performance  Computing  (HPC)  based  next  generation 
comprehensive  rotorcraft  analysis.  From  this  survey  of  the  key  technologies  in  other  application  areas  it  is 
possible  to  identify  the  critical  technology  gaps  that  stem  from  unique  rotorcraft  requirements. 


INTRODUCTION 

This  paper  presents  a  survey  of  computational  aeroe- 
lasticity  in  the  disciplines  of  fixed  wing  aircraft,  turboma¬ 
chinery,  and  rotary  wing  aircraft.  The  work  was  under¬ 
taken  by  the  U.S.  Army  Aeroflightdynamics  Directorate 
under  the  High  Performance  Computing  Institute  for  Ad¬ 
vanced  Rotorcraft  Modeling  and  Simulation  (HI- ARMS) 
and  NASA. 

The  survey  covers  the  emerging  methods  which  in¬ 
tegrate  Reynolds- averaged  Navier-Stokes  (RANS)  CFD, 
Finite  Element  Method  (FEM)  based  structural  mechan¬ 
ics,  and  high-fidelity  coupling  procedures  designed  to  sat¬ 
isfy  the  unique  requirements  of  each  discipline.  In  each 
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discipline,  the  key  aeromechanical  phenomena  which  de¬ 
termine  the  costs  and  risks  associated  with  design,  but 
remain  beyond  current  prediction  capabilities,  are  de¬ 
scribed.  The  current  status  of  High  Performance  Com¬ 
puting  (HPC)  based  high-fidelity  studies  on  the  key  phe¬ 
nomena  are  reviewed,  the  recent  advances  summarized, 
and  the  unresolved  challenges  highlighted. 

The  central  theme  of  this  paper  is  rotorcraft.  The 
objective  of  the  authors  —  together  with  a  companion  pa¬ 
per  on  requirements  [1]  —  is  to  lay  out  a  path  for  a  High 
Performance  Computing  (HPC)  based  truly  comprehen¬ 
sive  next  generation  rotorcraft  code;  comprehensive  in 
solutions  (performance,  loads,  stability,  vibration,  han¬ 
dling  qualities),  comprehensive  in  applications  (ground, 
hover,  steady  flight,  unsteady  maneuvers),  and  compre¬ 
hensive  in  scope  (arbitrary  geometries,  innovative  con¬ 
figurations).  The  intention  of  the  present  review  is  to 
identify  the  key  technologies  in  other  application  areas 
that  can  be  drawn  upon  to  this  end,  and  to  identify  the 
critical  technology  gaps  that  stem  from  unique  rotary 
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wing  requirements. 

The  paper  is  divided  into  five  sections.  The  first 
section  presents  a  description  of  the  state-of-the-art  in 
high-fidelity  fluid  structure  coupling.  This  is  followed  by 
three  sections,  one  each  on  the  status  of  computational 
aeroelasticity  in  fixed  wing  aircraft,  turbomachinery,  and 
rotary  wing  aircraft.  Each  section  is  subdivided  into  two 
parts;  the  first  part  is  on  structural  mechanics,  the  sec¬ 
ond  part  is  on  coupled  fluid-structure  applications.  The 
CFD  methods  applicable  to  each  discipline  are  not  re¬ 
viewed  in  this  paper.  The  last  section  summarizes  the 
different  CFD/CSD  coupling  nomenclatures  used  in  the 
three  disciplines. 

FLUID-STRUCTURE  COUPLING 

Definition  of  the  problem 

The  problem  of  fluid-structure  coupling  involves 
three  issues:  (1)  temporal  accuracy,  (2)  spatial  accuracy, 
and  (3)  interface  geometry  representation.  The  purpose 
of  the  present  discussion  is  to  clarify  their  meaning,  ex¬ 
plain  why  they  were  of  little  importance  in  rotorcraft  so 
far,  and  highlight  why  they  are  important  now  and  for 
the  future. 

At  the  PDE  level,  the  three  issues  are  related  to  the 
following  questions.  The  first  issue,  that  of  temporal  ac¬ 
curacy,  is  related  to  the  question  whether  two  systems 
of  coupled  PDEs  of  different  types  can  be  solved  sepa¬ 
rately,  one  after  the  other,  one  step  at  a  time,  exchang¬ 
ing  solutions  at  each  time  step.  This  is  the  method  of 
partitioned  formulations  that  begins  with  an  acceptance 
that  fluid  PDEs  of  convection-diffusion  type,  and  struc¬ 
tural  PDEs  of  elliptic  type  are  best  solved  separately  in 
their  own  domains  using  their  own  efficient  solvers.  The 
main  task  then  is  to  devise  a  method  of  solution  exchange 
which  renders  the  process  at  least  as  time  accurate  as  the 
worse  of  the  two  solvers.  The  second  and  the  third  issues, 
those  of  spacial  accuracy  and  interface  geometry  repre¬ 
sentation,  are  related  to  errors  introduced  during  solution 
exchange  across  domains.  Depending  on  these  errors  the 
method  of  solution  exchange  must  be  re-constructed  to 
achieve  an  intended  temporal  accuracy.  Note  that  the 
second  and  the  third  issues  are  independent  of  the  first, 
that  is,  they  arise  whether  or  not  a  partitioned  procedure 
is  adopted. 

Partitioned  vs.  fully  coupled  formulations 

The  practical  benefits  of  partitioned  formulations 
are  obvious  —  modularity  of  framework,  and  domain  sep¬ 
aration  for  refinements.  However,  an  appropriate  method 
of  solution  exchange  must  be  designed  to  ensure  time  ac¬ 
curacy.  Fully  coupled  formulations,  on  the  other  hand, 
are  strictly  time  accurate  without  this  additional  bur¬ 


den.  However,  such  formulations  require:  (1)  a  com¬ 
mon  time  integrator,  and  (2)  solution  of  an  algebraic  sys¬ 
tem,  at  each  Newton-like  step,  that  includes  both  fluid 
and  structural  stiffness.  The  second  requirement  is  not 
easy  to  meet  for  the  following  reasons.  Direct  solution 
is  not  an  option,  since  modern  CFD  (without  the  use 
of  reduced  order  modeling)  provides  far  too  many  de¬ 
grees  of  freedom  (DOFs),  typically  10-100M  for  rotor- 
craft.  Iterative  solution  is  the  natural  option  in  CFD. 
The  temporal  evolution  of  a  convection-diffusion  equa¬ 
tion  is  naturally  analogous  to  the  iterative  convergence 
of  an  algebraic  system.  Moreover,  they  are  easy  to  par¬ 
allelize.  On  the  other  hand,  iterative  solvers  are  tradi¬ 
tionally  not  preferred  for  structures.  The  convergence 
rates  of  iterative  solvers  depend  on  the  condition  num¬ 
ber  of  the  system  (for  symmetric  positive  definite  ma¬ 
trices  the  condition  number  is  the  ratio  of  the  largest  to 
the  smallest  eigenvalues,  or  natural  frequencies  squared). 
Typically,  aerospace  structures  pose  4-th  order  elastic¬ 
ity  problems  involving  bending-torsion-extension  of  thin 
beams,  or  bending-shear-extension  of  thin  plates,  and 
shells.  Condition  numbers  for  such  structures  range  from 
108  to  109  (the  number  for  a  typical  rotor  beam  model 
is  around  109).  In  comparison,  the  maximum  condition 
number  for  a  2-nd  order  fluid  problem  can  be  as  high  as 
106.  Pre-conditioners  that  are  well  suited  for  structures 
are  not  suitable  for  fluids.  For  example,  the  incomplete 
Cholesky  pre-conditioner  cannot  be  applied  to  fluids  as 
the  fluid  system  is  not  symmetric  positive  definite.  Simi¬ 
larly,  pre-conditioners  well  suited  for  fluids,  like  the  block 
Jacobi  (easily  parallelizable) ,  can  only  be  used  for  struc¬ 
tures  with  special  re-ordering  and  sparsity  fill-in  for  rea¬ 
sonable,  yet  problem  dependent,  convergence.  Moreover, 
such  procedures  demonstrate  poor  scalability.  Thus,  the 
challenge  of  a  unified  solution  procedure  in  the  presence 
of  CFD  is  real  [2]. 

The  difficulty  is  bypassed  for  structures  that  are  less 
ill-conditioned.  One  example  is  turbomachinery  struc¬ 
tures.  These  are  modeled  using  solid  elements,  gov¬ 
erned  by  2-nd  order  elasticity.  Gottfried  and  Fleeter 
[3]  have  used  full  coupling  for  turbomachinery  aeroelas¬ 
ticity.  Their  analyses,  TAM-ALE3D,  a  refined  version 
of  the  original  ALE3D  code  developed  in  the  Lawrence 
Livermore  National  Laboratories,  is  a  fully  coupled  finite 
element  analysis  that  has  been  applied  to  turbomachin¬ 
ery  flutter  calculations.  A  second  example  is  biological 
structures.  Here,  structures  are  either  membrane- like,  or 
thick,  when  shell-like.  Bathe  and  Zhang  [4]  have  applied 
full  coupling  in  biomedical  fluid  flows,  as  well  as  general 
internal  flows  in  mechanical  components.  Their  analysis 
is  used  extensively  for  biomedical  and  mechanical  appli¬ 
cations,  as  part  of  the  ADINA  commercial  software  code. 
An  option  for  partitioned  formulation  is  also  provided 
(the  authors  denote  them  as  direct  and  iterative  proce- 
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dures).  The  foundations  of  all  fully  coupled  approaches 
can  be  traced  back  to  the  Arbitrary  Lagrangian  Eulerian 
(ALE)  formulation  of  conservations  laws,  first  introduced 
by  Hirt  et  al.  [5]  in  1974  for  fluid  flows  at  all  speeds.  It 
was  applied  by  Donea  [6]  in  1983  for  fluid-structure  in¬ 
teraction  problems.  Bendikson  [7,  8]  first  applied  this 
technique  in  the  early  1990s  for  flutter  calculations  and 
demonstrated  minimal  errors  in  energy  transfer  between 
fluids  and  structures.  Fully  coupled  formulations  are  also 
termed  monolithic  formulations. 

In  summary,  for  aerospace  structures,  partitioned 
formulations  provide  fundamental  advantages  over  fully 
coupled  ones,  in  addition  to  their  obvious  practical  ben¬ 
efits. 

Current  coupling  practices  in  rotorcraft 

Classical  comprehensive  analyses,  which  use  lifting¬ 
line  or  lifting-surface  models,  involve  at  the  most  1-10K 
DOFs.  The  formulation  is  fully  coupled,  and  is  solved 
using  direct  methods.  Thus  the  first  issue  of  temporal 
accuracy  is  satisfied.  To  re-iterate,  the  need  for  a  par¬ 
titioned  formulation  for  fluids  and  structures  is  felt  only 
under  the  following  two  circumstances.  First,  when  mil¬ 
lions  of  fluid  DOFs  call  for  an  iterative  solver  whereas  the 
structural  DOFs,  being  ill-conditioned,  call  for  a  direct 
solver,  or  at  least  a  different  pre-conditioner.  Second, 
when  these  conflicting  requirements  of  fluids  and  struc¬ 
tures  are  very  easily  met  in  separate  domains.  The  sec¬ 
ond  issue  of  spatial  accuracy  is  easily  satisfied  by  such  for¬ 
mulations  as  the  structural  shape  functions  are  available 
to  the  fluid  terms.  The  third  issue  of  interface  geometry 
representation  has  also  been  simple,  so  far,  in  rotorcraft. 
This  is  because  the  beam  theory,  or  any  single  compo¬ 
nent  structural  model,  regardless  of  1-D,  2-D,  or  3-D, 
provides  a  simple  yet  rigorous  definition  of  the  surface. 
The  complexity  arises  in  multi-component  structures  like 
the  fuselage. 

Dynamics  researchers  have  carried  out  coupled 
rotor- fuselage  analysis,  but  at  a  time  when  CFD  air¬ 
loads  were  beyond  state  of  the  art.  With  the  emergence 
of  this  capability,  there  is  a  requirement  to  address  the 
issue  of  geometry  representation.  A  detailed  dynamic 
model  is  not  necessarily  the  best  suited  for  fluid-structure 
coupling.  The  surface  geometry  representation  is  more 
important  than  the  internal  load  paths.  Thus,  the  re¬ 
quirement  is  to  have,  at  the  least,  a  shell  model  that 
re-produces  the  low  frequency  modes  (up  to  40  Hz),  and 
at  best,  a  detailed  model  that  includes  the  outer  skin. 
Fixed-wing  researchers  have  already  accomplished  these 
tasks,  and  there  is  a  volume  of  literature  that  can  be 
drawn  upon. 

Rotorcraft  CFD/CSD  coupled  analyses  have  used 
partitioned  procedures.  The  focus  so  far  has  been  on  the 


rotor.  For  an  isolated  rotor,  the  wetted  surface  is  rigor¬ 
ously  defined.  The  shape  functions  are  available.  Thus 
the  issues  of  spatial  accuracy  and  surface  representation 
are  easily  dealt  with.  The  issue  of  temporal  accuracy  is 
enforced  concurrently  with  Newton-Raphson  type  itera¬ 
tions  for  determining  the  trim  angles.  For  a  time  march¬ 
ing  solution,  the  procedure  is  conceptually  simpler  as  the 
trim  angles  are  left  undetermined,  but  requires  the  en¬ 
forcement  of  temporal  accuracy  at  every  time  step.  The 
problem  is  then  the  same  as  that  faced  by  fixed-wing  re¬ 
searchers,  where  it  has  been  of  great  importance  due  to 
the  emphasis  on  flutter.  Small  errors  in  coupling  result  in 
erroneous  energy  exchange  and  affect  flutter  boundaries 
adversely.  The  emphasis  in  rotors  has  been  on  loads. 
Errors  in  fluid-structure  coupling  are  less  visible.  The 
problem  is  further  compounded  in  fixed-wing  because  of 
its  3-D  structure  and  multiple  sub-structures.  Without  a 
detailed  3-D  structure,  this  difficulty  has  not  been  faced 
yet  in  rotorcraft. 

Temporal  accuracy 

The  temporal  accuracy  of  partitioned  formulations 
is  the  main  focus  of  high-fidelity  fluid-structure  coupling. 
Physically,  it  means  that  the  airloads  and  structural 
loads  at  a  given  time  step  are  consistent  with  one  an¬ 
other.  Numerically,  it  means  that  the  temporal  error  is 
driven  down  to  the  level  of  the  worse  solver. 

One  effective  approach  is  the  use  of  sub-iterations 
between  every  consecutive  pair  of  time  steps.  The 
method  was  first  applied  in  aeroelastic  computations  by 
Strganac  and  Mook  [9] ,  followed  by  Weeratunga  and  Pra- 
mono  [10],  and  more  recently  by  Melville  and  his  co¬ 
researchers  [11,  12].  The  method  is  time  intensive  for 
RANS.  Potential  innovations  may  involve  the  coupling 
of  structural  dynamics  with  the  fluid  sub- iterations. 

Early  work  in  fixed-wing  fluid  structure  coupling  in¬ 
volved  sub-iteration  free,  straight-forward  time  integra¬ 
tion.  See  for  example,  Edwards  [13],  Bennet  [14],  and 
Guruswamy  et  al.  [15,  16].  In  the  serial  approach,  one 
solver  waited  for  the  completion  of  the  other,  before  ex¬ 
changing  solutions  and  advancing  to  the  next  time  step. 
In  the  parallel  approach,  both  solvers  advanced  simulta¬ 
neously,  followed  by  solution  exchange  before  advancing 
to  the  next  step.  The  latter  was  proposed  originally  by 
Weeratunga  and  Pramono  [10],  and  refined  subsequently 
by  Far  hat  and  Lesoinne  [17].  The  original  work  involved 
accompanying  sub-iterations  due  to  the  relative  instabil¬ 
ity  of  the  scheme.  The  subsequent  refinement  improved 
stability  without  sub- iterations,  at  the  cost  of  exchang¬ 
ing  solutions  twice  at  each  time  interval.  Since  then,  a 
significant  amount  of  research  has  been  conducted  on  de¬ 
vising  sub- iteration  free  methods  that  retain,  at  the  least, 
second-order  accuracy.  Far  hat  and  his  co-researchers  [18] 
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have  described  the  formal  design  of  such  methods. 

Sub-iteration  free  approaches  are  predictor-corrector 
schemes,  tailored  to  the  individual  time  integrators.  A 
generalized  predictor-corrector  approach  is  illustrated 
below,  adopted  from  the  last  reference.  The  mesh  defor¬ 
mation  is  denoted  by  x,  the  fluid  variables  are  denoted  by 
q,  and  the  structural  deformations  by  u.  The  time  step 
is  denoted  by  n.  The  fluid-structure  interface  boundary 
is  denoted  by  T.  For  example  uy  denotes  the  structural 
deformations  at  the  interface,  a  subset  of  u.  The  symbol 
< —  denotes  a  time  integrator  and  shows  dependencies 
on  the  latest  time  steps.  A  single  time  step  advancement 
is  given  as  follows. 

1.  Predict  deflections  for  the  next  time  step 


5.  Update  structural  solver 

un+1  < — un,  F™+1  (including  up*  ) 

The  predictor  is  on  displacement  (step  1).  The  corrector 
is  on  forcing  (step  4).  Methods  of  the  above  type  are 
classified  as  loose  coupling  in  the  fixed  wing  community. 
If  the  above  procedure  is  carried  out  multiple  times  at  the 
same  time  step  (i.e.  sub- iterations),  then  the  following 
can  be  enforced 

u™+lP  =  u ™+1 

which  implies  that  the  mesh  deformation  follows 
4+1  =r?  +  Ur  «+1  -  u£) 


u 


n+ 1 

r 


p 


2.  Update  mesh  points  x.  Mesh  boundary  points  de¬ 
noted  by  xy •  Mesh  internal  points  denoted  by  xq. 
The  boundary  mesh  points  can  be  updated  as 


x ™+i 


p 


n+P 


For  example, 

xn+lF  =  xf  +  Ur  (uplP  -  <P) 


The  above  equation  simply  represents  the  discretized 
deflection  and  velocity  compatibilities  on  the  sur¬ 
face.  Without  re-griding  in  time,  the  deflection  com¬ 
patibility  is  given  by  x  =  Up^  where  Up  is  simply 
the  transformation  that  connects  the  CFD  surface 
grid  to  the  structural  grid.  For  perfectly  matched 
meshes  Ur  is  equal  to  the  identity  matrix.  The  ve¬ 
locity  compatibility  is  x  =  Uph.  If  the  structure 
reaches  out  to  the  wetted  surface,  and  the  shape 
functions  are  known,  Up  =  H,  where  H  are  the 
shape  functions.  With  re-griding,  the  mesh  connec¬ 
tivity  changes  with  time,  i.e.  the  transformation  Up 
is  now  a  function  of  time,  and  the  discretized  veloc¬ 
ity  compatibility  is  expressed  as  x  =  Up  ft  where  Up 
denotes  a  selected  combination  of  Up  over  previous 
time  steps,  e.g.  the  mean  of  Up  and  Up-1. 


The  internal  points  are  updated  based  on  the  bound¬ 
ary  points 


rn+U 


rU  n+1 


3.  Update  fluid  solver  using  updated  mesh 

qn+1  qn,xn+lP 


rather  than 

£p+1  =  x7^  +  Up  ^p+1  —  Uy  ^ 

and  hence  enforces  the  velocity  compatibility  x  =  Chip 
strictly.  The  forcing  and  structural  response  (step  4  and 
5)  are  also  consistent,  as  Xp+1  is  replaced  with  Xp+1. 
The  method  with  sub-iterations  is  classified  as  close  cou¬ 
pling,  tight  coupling ,  or  strong  coupling  in  the  fixed  wing 
literature.  The  force  predictor  in  step  4  can  be  chosen  in 
various  ways.  The  first  form  is  also  called  serial  stagger¬ 
ing.  The  second  form  is  called  parallel  staggering,  simply 
because  the  fluid  and  structural  updates  (step  3  and  5) 
can  be  performed  independently  of  one  another. 

The  methods  that  enforce  second-order  accuracy, 
without  sub- iterations,  rely  on  the  formal  selection  of: 
(1)  the  integrators,  i.e.  the  arrows  b — (2)  the  initial 
predictor,  i.e.  step  1,  and  (3)  the  form  of  the  force  cor¬ 
rector,  i.e.  step  4.  Two  illustrative  examples  are  provide 
in  reference  [18]  based  on  a  second-order  fluid  integra¬ 
tor,  and  a  second-order  (Newmark  type)  structural  in¬ 
tegrator.  Unless  formally  selected,  the  sub-iteration  free 
methods  provide  only  first-order  accuracy. 

Lastly,  we  note  that  for  non-CFD  Lagrangian  aero¬ 
dynamic  models,  the  time  iteration  procedure  is  simpler 
as  there  are  no  grid  motions.  The  need  to  reconcile  the 
Geometric  Conservation  Law  and  structural  motion  up¬ 
date  does  not  arise.  The  structural  forcing  Fs  is  related 
directly  to  the  deformations  u.  The  serial  and  parallel 
staggered  schemes  along  with  sub-iterations  can  be  im¬ 
plemented  as  follows. 

1.  Serial  staggered  with  sub- iterations 
Predict  un+lP 


4.  Calculate  structural  forcing 


pn+l 


=  Fs(qn+\xplP) 

or  Fs  (qn,XY  V..  or  other  forms. 


Then,  perform  successively 

pn+i  4 —  FJ\itra+lP;  in  step  1 

un+ 1  < —  rd\Fsn+1;  in  step  2 

Iterate  until  convergence,  i.e.  unJrlP  ~  unJrl 
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2.  Parallel  staggered  with  sub- iterations 
Predict  un+lP 

Then,  perform  in  a  single  step 

Fn+1  < pn^un+ 1P. 

un+l  i -  un^Fn. 

Exchange  updates,  i.e.  update  un+lP  with  un+1  in 
the  airloads  integrator,  and  update  F'1  with  F"+l 
in  the  structures  integrator,  and  iterate  until  con¬ 
vergence. 

The  parallel  staggered  method  is  advantageous  when 
the  fluid  and  the  structure  are  run  on  separate  processors. 
The  parallel  method  has  inferior  stability  and  requires 
requires  lower  time  steps. 

Spatial  accuracy 

Assume  that  the  structural  model  includes  a  rig¬ 
orous  interface  geometry  representation.  However,  the 
CFD  and  the  CSD  discretizations  will  not,  in  general, 
match  at  the  interface.  The  deformation  and  struc¬ 
tural  forcing  must  be  mapped  correctly  across  the  non¬ 
matching  interface.  Generic  methods  which  are  ‘exact’ 
regardless  of  the  blade  shape  are  critical  for  the  evalu¬ 
ation  of  advanced  geometry  blades.  Two  such  methods 
are  formulated  below. 

Deformation  mapping  deals  with  the  correct  evalua¬ 
tion  of  Ur  in  step  2  (see  previous  section).  Loads  transfer 
deals  with  the  correct  evaluation  of  structural  forcing  Fs 
in  step  4.  The  latter  implies  an  exact  evaluation  of  the 
virtual  work  term.  This  is  a  necessary  energy  conserva¬ 
tion  requirement.  In  addition,  it  is  desirable,  though  not 
necessary,  that  the  total  integrated  forcing  be  preserved 
during  the  mapping. 

Deformation  mapping  is  straight-forward  when  the 
shape  functions  are  available.  Loads  transfer  can  be  ac¬ 
complished  via:  (1)  integrated  force  coupling,  and  (2)  di¬ 
rect  traction  coupling.  In  the  integrated  force  coupling, 
the  virtual  work  is  calculated  based  on  integrated  fluid 
stresses  (pressure  and  skin  friction)  over  a  surface  patch. 
In  surface  traction  coupling,  the  virtual  work  is  calcu¬ 
lated  based  on  the  pressure  and  shear  distributions  di¬ 
rectly.  The  first  method  preserves  total  forces,  but  does 
not  guarantee  energy  conservation  for  a  finite  mesh.  The 
second  method  is  strictly  energy  conservative,  but  does 
not  guarantee  preservation  of  total  forces.  Thus,  they 
are  complimentary  to  each  other.  However,  as  long  as  in¬ 
terpolations  and  integrations  are  performed  consistently 
within  each  domain,  both  satisfy  and  conservation  and 
preservation  in  the  limit  of  mesh  refinement. 

Both  methodologies  can  be  formulated  for  rotor 
blades.  Simple  illustrations  are  given  below. 


F 

Vx 


Figure  1:  A  concentrated  force  on  a  rotor  blade  ob¬ 
tained  by  integrating  fluid  stresses  over  a  surface 
patch  of  area  A  A 

The  integrated  force  coupling  is  formulated  as  fol¬ 
lows.  Consider  a  concentrated  force  F  acting  at  a  point 
P  on  the  deformed  blade,  Fig.  1.  F  is  integrated  trac¬ 
tion  over  an  arbitrary  patch  of  area  A  A.  The  virtual 
work  done  by  the  force  is  simply 

SW  =  F  •  SrP  (1) 

where  rp  is  the  position  vector  of  the  point  P,  and  Srp  is 
its  virtual  displacement.  The  position  can  be  expressed 
as 

rP  =  rTe°  (2) 

r  are  the  scalar  components  and  e°  =  [ijk]T  is  a  set  of 
unit  vectors.  Let  u  be  the  states  of  blade  deformation. 
The  virtual  displacement  then  becomes 

Srp  =  (DSu)T  e°  (3) 

where  D  is  the  derivative  matrix  of  the  scalar  components 
r  with  respect  to  the  states  u.  Henceforth,  matrices  are 
denoted  in  bold  capital  letters.  In  beam  theory,  the  states 
are  typically  three  deflections  and  three  rotations  u  = 
[ui ,  U2 ,  ^3,  #i,  #2?  #3]-  F)  is  then  the  (3  x  6)  derivative 
matrix.  Expressing  the  force  in  the  same  basis 

F  =  FTe°  (4) 

gives 

SW  =  FtD5u  (5) 

If  q  are  the  N  generalized  nodal  displacements  of  a  finite 
element  containing  the  point  P,  and  H  is  the  (6  x  N) 
elemental  shape  function  matrix,  it  follows  from  Su  = 
H  Sq 

SW  =  FTDH<5q  =  QTSq  (6) 

The  generalized  nodal  force  is  then  given  by 

Q  =  UtBtF  (7) 

The  term  DT  transmits  the  airloads  in  3-D  space  to  the 
1-D  beam  structure.  The  matrix  D  varies  with  the  choice 
of  beam  theory.  Note  that,  from  eqn.  3,  the  virtual  dis¬ 
placement  components  are 

Srp  =  D  Su  =  DH  Sq  (8) 
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Equations  7  and  8  highlight  the  well-known  relation 
that  the  generalized  structural  forcing  vector  relate  to 
the  aerodynamic  forcing  via  the  transpose  of  the  relation 
that  connects  the  aerodynamic  deflections  to  structural 
deflections.  Here,  DH  can  be  interpreted  as  the  equiva¬ 
lent  elemental  shape  functions  in  3-D  space  for  the  cor¬ 
responding  beam  theory. 

Note  that  the  force  method  ensures  that  the  total 
integrated  forces  (and  moments)  remain  the  same  when 
transferred  from  the  fluid  to  the  structural  domain.  How¬ 
ever,  because  the  point  of  application  of  F  within  each 
surface  patch  is  arbitrary,  the  method  is  energy  conser¬ 
vative  only  in  the  limit  of  mesh  refinement.  Strict  con¬ 
servation  can  be  enforced  using  traction  coupling. 


Figure  2:  Fluid  pressure  and  surface  shear  over  a 
rotor  blade  differential  area  dA 

The  direct  traction  coupling  is  formulated  as  follows. 
Consider  a  differential  area  dA  within  the  surface  patch 
A  A  of  magnitude  dA  and  unit  normal  n,  Fig.  2. 

dA  =  ndA  =  nTe°dA  (9) 

The  differential  force  generated  by  the  pressure  is 

—pdA  =  —pnTe°dA 

The  differential  force  generated  by  the  fluid  stress  tensor 
along  the  direction  of  dA  is 

<tf  •  dA  =  (j f  •  ndA  =  (<tf n)Te°dA 


In  the  force  coupling  case,  the  following  integration  was 
performed  in  the  fluid  domain. 


F  = 


-j  - 


p  dA  +  <tf  •  dA 


/  pnTdA+  {ctfu)1 
Jaa  Jaa 


dA 


e°  =  FTe° 


(10) 


The  integrated  force  was  then  transferred  to  the  struc¬ 
tural  domain.  H  and  D,  required  for  this  calculation, 
were  available  in  the  structural  domain.  Now  consider 


the  virtual  work  calculation  using  fluid  traction. 


SW  =  j  (—pdA  +  crp  •  dA)  •  Srp 
Jaa  '  ' 

=  /  [—  pnT  +  ((7Fn)T]  e°  •  (D5u)T  e°dA  (11) 

Jaa 

=  /  [—  pnT  +  (dpn)7]  DH£gcL4  =  QT5q 

Jaa 

The  generalized  nodal  force  then  becomes 


Q  = 


j 

Jaa 


— HTDTnp  +  Hj  Dj  uf  n  \  dA 


Tt\T, 


(12) 


Note  that  the  integration  involves  variables  from  both 
domains.  If  performed  in  the  fluid  domain,  exact  values 
of  D  and  H  must  be  received  from  the  structural  do¬ 
main  at  the  fluid  Gauss  points.  Similarly  if  performed  in 
the  structural  domain,  exact  values  of  p  and  ap  must  be 
received  from  the  fluid  domain  at  the  structural  Gauss 
points.  The  value  of  n  is  different  in  both  (except  in  the 
ideal  case  when  the  meshes  and  spatial  orders  match  ex¬ 
actly).  Unlike  integrated  force  coupling,  direct  traction 
coupling  ensures  the  exact  calculation  of  virtual  work. 

There  has  been  significant  contributions  by  vari¬ 
ous  researchers  to  address  the  issue  of  conservation  and 
preservation,  within  the  context  of  temporal  accuracy. 
See  for  example,  Maman  and  Farhat  [19],  Cebral  and 
Lohner  [20],  Farhat  et  al.  [21],  Slone  et  al  [22],  and  Mich- 
ler  et  al.  [23].  The  key  conclusion  is  that  exact  conser¬ 
vation  and  preservation  can  be  ensured  in  the  limit  of 
mesh  refinement  only  when  interpolations  of  all  the  vari¬ 
ables  are  performed  using  schemes  consistent  with  their 
domain.  For  example,  p  and  ap  can  be  interpolated  only 
in  the  fluid  domain  and  then  transferred.  Similarly,  de¬ 
formations  can  be  interpolated  only  in  the  structural  do¬ 
main  and  then  transferred. 

Neither  of  these  methods  have  been  applied  to  rotor- 
craft  CFD/CSD  so  far.  In  the  current  applications,  the 
rotor  blade  is  excited  by  force  and  moment  distributions 
(per  span)  along  its  local  sectional  elastic  axis.  The  force 
and  moment  excitations  are  obtained  by  integrating  the 
fluid  stresses  along  chord- wise  strips.  Chord- wise  strips 
are  an  easy  and  natural  choice  for  structured  grids.  The 
excitations,  if  calculated  at  the  structural  Gauss  points 
along  the  elastic  axis  (using  interpolation  in  the  fluid  do¬ 
main),  or  imposed  as  concentrated  forces,  one  for  each 
strip  using  the  shape  functions,  would  satisfy  conserva¬ 
tion  and  preservation  requirements  in  the  limit  of  mesh 
refinement. 

The  exact  generic  methods  provide  significant  ad¬ 
vantages  over  the  current  methods  of  sectional  airloads. 
First,  the  method  of  sectional  airloads  is  arbitrary  for 
advanced  geometry  blades,  a  simple  example  of  which  is 
the  BERP  tip.  For  such  geometries  it  is  necessary  to 
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Figure  3:  FEM  model  (MSC  NASTRAN)  of  a  SH- 
60  Sea  Hawk  fuselage  developed  by  Sikorsky  and 
used  in  UMARC  rotor-fuselage  coupled  dynamic 
analysis;  Yang  et  al.  2004 


Figure  4:  FEM  model  of  a  F-16  fighter  used  in  a 
doublet-lattice  aerodynamic  model  based  typical 
fixed-wing  flutter  analysis;  Denegri  2005 

use  the  exact  generic  methods  for  high-fidelity  coupling. 
Second,  the  generic  methods  are  equally  applicable  for 
both  structured  and  unstructured  grids,  the  latter  being 
increasingly  applied  today  for  near-blade  flows.  Third, 
the  methods  are  well-suited  for  the  long  term  CFD  goal 
of  near  blade  adaptive  refinement.  Fourth,  the  generic 
methods  are  equally  applicable  for  the  rotor,  fuselage, 
and  any  structure  in  general,  as  long  as  the  interface  ge¬ 
ometry  representation  is  clearly  defined. 

Interface  geometry  representation 

The  issue  of  interface  geometry  representation  in¬ 
volves  the  rigorous  description  of  the  surface  based  on  the 
underlying  structural  model.  The  method  of  description 
differs  based  on:  (1)  whether  the  structural  shape  func- 


Figure  5:  Detailed  FEM  model  of  a  F-16  fighter 
with  0.17  M  DOFs  for  RANS  CFD/CSD  based 
aeroelastic  response  prediction;  Farhat  et  al. 
2003 

tions  are  available,  and  (2)  even  if  the  shape  functions  are 
available  whether  the  structural  model  reaches  out  to  the 
wetted  surface.  The  problem  posed  by  the  unavailabil¬ 
ity  of  shape  functions  is  a  practical  complication  that  is 
ideally  avoided.  The  problem  posed  by  the  structure  not 
reaching  out  to  the  wetted  surface  is  more  fundamen¬ 
tal.  For  a  single  component  structure,  e.g.  a  rotor  blade 
beam  model,  the  necessary  interpolation  and  extrapola¬ 
tions  are  rigorously  defined  by  the  underlying  theory  (e.g. 
note  the  term  D  in  eqn  7  in  addition  to  H).  Complica¬ 
tions  arise  for  multiple  intersecting  sub-structures,  as  is 
the  case  for  airframes. 

For  illustration,  figures  3  and  4  show  state-of-the- 
art  dynamic  models  of  rotary  and  fixed  wing  airframes. 
The  details  of  the  work  are  given  later,  in  the  respec¬ 
tive  survey  sections.  The  emphasis  in  dynamic  models 
is  on  the  key  load  bearing  components  and  internal  load 
paths.  The  emphasis  in  fluid-structure  coupling  is  on  the 
external  shell.  For  such  purposes  models  which  include 
the  outer  shell,  and  in  addition,  the  internal  load  paths, 
are  most  appropriate.  For  example,  Fig.  5,  shows  a  more 
detailed  model  of  the  same  F-16,  including  the  outer  shell 
(details  of  the  work  is  cited  later). 

Detailed  modeling  of  the  fuselage  for  fluid-structure 
interaction  has  not  been  a  key  focus  so  far  in  rotorcraft 
for  the  following  reason.  The  dominant  dynamics  of  a 
helicopter  fuselage,  unlike  fixed-wing,  is  vibration,  which 
occurs  at  high  frequencies  (usually  pN^/ rev).  The  domi¬ 
nant  source  of  this  vibration  is  not  the  fuselage  flow  field, 
but  the  shaft  transmitted  forcing  from  the  rotor.  The 
forcing  from  the  rotor  primarily  stems  from  the  rotor 
flow  field,  the  interactional  effect  of  the  fuselage  on  the 
rotor  flow  field  occurs,  nominally,  at  lower  frequencies 
(1/rev). 

There  are  important  conditions  where  the  above  ar¬ 
gument  is  invalid.  First,  at  low  speed  flight  (around 
/i  =  0.15)  direct  impingement  of  the  rotor  tip  vortices 
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(a)  1-D  line  fuselage 


(b)  3-D  detailed  fuselage 

Figure  6:  Measured  and  NASTRAN  calculated 
fuselage  natural  frequencies  of  the  AH-1G  heli¬ 
copter;  Yeo  and  Chopra  2002 

on  the  tail  boom  is  a  key  mechanism  of  vibration  (at 
Nb/ rev).  RANS  was  unable  to  resolve  this  flow  field  thus 
far,  but  is  beginning  to  do  so  now.  Second,  for  advanced 
configurations  with  low  shaft  clearance  (for  low  drag)  di¬ 
rect  wake  impingement  is  a  significant  source  of  vibration 
for  a  wide  range  of  flight  speeds.  There  has  been  limited 
resolution  of  the  interactional  flow  field  around  the  fuse¬ 
lage  thus  far.  RANS  based  CFD,  and  alternative  wake 
CFD  models  like  the  Vorticity  Transport  Model  [24,  25] , 
are  beginning  to  do  so.  The  emergence  of  interactional 
CFD  opens  opportunity  to  calculate  vibration  and  buffet 
loads  in  response  to  the  interactional  flow  field. 

Reproducing  the  structural  frequencies  of  the  fuse¬ 
lage  (need  at  least  up  to  30  Hz)  is  in  itself  difficult  in 
helicopters  (more  than  20%  error  without  detail  model¬ 
ing).  Difficult  components  he  in  critical  load  paths  (main 
rotor  hub/pylon/transmission  case).  Dynamicists  often 


(a)  1-D  line  fuselage 


(b)  3-D  detailed  fuselage 

Figure  7:  FEM  line  and  3-D  fuselage  NAS¬ 
TRAN  models  of  the  AH-1G  helicopter  from  the 
DAMVIBS  program  used  for  rotor-fuselage  cou¬ 
pled  dynamic  analysis;  Yeo  and  Chopra  2002 


use  line  fuselages  to  re-produce  measured  frequencies. 
For  example,  the  calculated  frequencies  of  an  AH-1G 
(two-bladed,  teetering)  fuselage  for  two  different  types 
of  models  are  compared  with  measurements  in  Fig.  6. 
The  corresponding  models  are  shown  in  Fig.  7.  The  sim¬ 
pler  model  is  designed  to  generate  similar  frequencies  as 
the  more  detailed  model,  and  is  often  preferred.  The 
fluid- structure  coupling  methodology,  however,  is  more 
involved  for  the  simpler  model  compared  to  the  detailed 
model.  The  need  for  extrapolation,  a  key  source  of  error, 
is  minimized  in  the  case  of  detailed  models  that  include 
the  external  shell  structure.  Models  with  both  internal 
and  external  details  of  the  structure  are  used  regularly  by 
crash  researchers,  for  very  short-time  high  impact  tran¬ 
sients  (less  than  0.1  sec).  An  example  of  a  crash  test 
simulation  is  shown  in  Fig.  8,  taken  from  Ref.  [26]. 

Constructing  a  generalized  modular  interface  which 
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(a)  Detailed  shell  fuselage 


Acceleration,  g 


(b)  Measured  and  predicted  acceler¬ 
ations  at  a  front  left  fuselage  station 

Figure  8:  Detailed  FEM  composite  fuselage  (MSC 
Dytran)  as  used  in  high-fidelity  crash  models; 
Fasanella  et  al.  2007  (courtesy  Karen  Jackson, 
NASA  Langley) 


is  adaptable  to  alternative  levels  of  underlying  struc¬ 
tural  fidelity  requires  innovative  interpolation  and  ex¬ 
trapolation  methods  for  interface  displacement  mapping. 
A  large  volume  of  literature  has  been  devoted  to  such 
methods  for  fluid-structure  coupling.  Broadly,  they  are 
divided  into  two  categories.  First,  the  surface  tracking 
methods.  These  rely  on  the  underlying  shape  functions 
while  trying  to  mitigate  the  complexities  associated  with 
multi-component  sub- structures.  Second,  the  surface  fit¬ 
ting  methods.  These  are  focused  more  on  independent 
(from  fluid  and  structural  solvers)  formulations  without 
the  need  for  shape  functions.  In  fixed  wing  applications, 
the  current  practice  is  to  use  a  combination  of  these  in¬ 
terpolation  methods  to  map  deformations  from  CSD  to 


CFD.  Mapping  airloads  from  CFD  to  CSD  then  simply 
amounts  to  the  correct  evaluation  of  virtual  work  based 
on  the  deformation  mapping  (compare  eqn.  8  with  eqns.  7 
and  12). 

A  displacement  mapping  must:  (1)  recover  large  de¬ 
flections  and  rotations  of  the  underlying  structure  from 
the  surface  displacement,  and  (2)  ensure  displacement  as¬ 
sociation  with  time,  i.e.  attachment  points  must  remain 
attached  throughout  the  simulation.  Arbitrary  proce¬ 
dures  like  a  nearest  neighbor  mapping,  although  simple 
and  fast,  do  not  meet  these  requirements.  Moreover,  it 
prevents  the  correct  evaluation  of  the  virtual  work  term. 
Procedures  which  rely  on  underlying  structural  shape 
functions  for  multi-component  blending  and  extrapola¬ 
tion  meet  these  requirements.  They  are  called  surface 
tracking  methods  as  mentioned  earlier.  When  the  shape 
functions  are  not  available,  the  most  general  procedure 
is  to  use  radial  basis  functions.  These  are  also  called 
surface  fitting  methods.  The  well-known  surface  fit¬ 
ting  methods  include  the  Multiquadric-Biharmonic  (MQ) 
method,  the  method  of  surface  splines  commonly  called 
the  Infinite  Plate  Spline  (IPS),  the  method  of  Thin 
Plate  Spline  (TPS),  and  Non-Uniform  Rational  B-Spline 
(NURBS).  These  methods  have  been  systematically  stud¬ 
ied  by  Hounjet  and  Meijer  [27]  in  1994,  and  Smith  et 
al.  [28]  in  2000  (see  references  therein  for  details  of  the 
above  methods).  The  study  by  Smith  et  al.  found  TPS 
to  be  most  suitable  for  interpolation  and  extrapolation 
of  generalized  structural  geometries.  These  radial  ba¬ 
sis  functions  are  defined  on  the  whole  domain  of  a  sub¬ 
structure.  Complications  arise  for  multiple  intersecting 
sub-structures  in  3-D  space. 

Fitting  methods  applicable  in  3-D,  which  do  not  use 
radial  basis  functions,  have  also  been  implemented,  no¬ 
table  among  which  are  the  Constant  Volume  Tetrahe¬ 
dral  (CVT)  method  by  Goura  et  al.  [29,  30],  and  the 
Inverse  Boundary  Element  (IBE)  method  by  Chen  and 
Jadic  [31].  Even  though  applicable  in  3-D,  these  methods 
are  less  suitable  for  large  deformation  nonlinear  prob¬ 
lems.  The  first  method  provides  a  nonlinear  mapping. 
Linearization  destroys  the  exactness  of  rigid  body  rota¬ 
tions.  The  second  method  requires  iterations  to  conform 
to  the  deformed  geometry,  unless  the  deformations  are 
assumed  to  he  within  the  range  of  linear  elasticity. 

An  emerging  technique  is  the  use  of  3-D  radial  basis 
functions  with  compact  support.  These  are  defined  lo¬ 
cally,  and  hence  are  easily  applicable  for  surfaces  with 
high  3-D  curvature.  The  method  is  based  on  multi¬ 
variate  scattered  data  interpolation.  Originally  proposed 
by  Wu  [32]  and  Wendland  [33] ,  and  subsequently  applied 
by  Beckert  and  Wendland  [34]  for  fluid-structure  inter¬ 
action  problems,  the  method  has  found  increasing  appli¬ 
cation  in  recent  large-scale  fixed-wing  CFD/CSD  calcu¬ 
lations  [35,  36]. 
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FIXED  WING  AIRCRAFT 

A  comprehensive  survey  of  the  state-of-the-art  in 
fixed-wing  aeroelasticity  can  be  found  in  Livne  [37]  in 
2003.  Schuster  et  al.  [38],  in  the  same  year,  have  reviewed 
the  state-of-the-art  in  computational  aeroelasticity.  The 
term  was  used  broadly,  covering  unsteady  aerodynam¬ 
ics  to  CFD  and  beam  models  to  detailed  FEM.  In  2004, 
Kamakoti  and  Shyy  [39]  described  the  RANS/FEM  cou¬ 
pling  procedures  used  for  fixed  wing  applications.  The 
AGARD  445.6  wing  as  taken  as  the  baseline  configura¬ 
tion  for  illustration.  More  recent  developments  in  com¬ 
putational  aeroelasticity  are  highlighted  by  Bartels  and 
Sayma  [40]  in  2007.  In  the  same  year,  a  report  detailing 
European  research  towards  the  development,  evaluation, 
and  transition  of  high-fidelity  aeroelastic  tools  for  full 
aircraft  has  also  been  published  [41]. 

The  following  sections  review  the  status  of  compu¬ 
tational  aeroelasticity  in  fixed  wing  applications.  The 
first  section  briefly  describes  the  status  of  detailed  CSD. 
The  second  section  summarizes  the  important  aeroelas¬ 
tic  phenomena  and  the  status  of  high-fidelity  CFD/CSD 
methods  in  predicting  them. 

Structural  Modeling 

First,  the  meaning  of  ‘detailed  structures’  must  be 
clarified.  Shown  in  Fig.  9  is  a  typical  Boeing  737  wing  sec¬ 
tion  with  a  simple  main/aft  double  flap  system  with  a  sin¬ 
gle  slotted  thrust  gate,  in  its  retracted  position  [42].  Only 
the  trailing-edge  is  shown,  the  three  position  Kruger  slats 
in  the  leading  edge  are  not  shown.  Typical  structures 
such  as  this,  incorporate  kinematic  constraints  to  rigid 
body  modes  to  flexible  deformations.  For  aeroelastic  and 
aeroservoelastic  studies,  a  distinction  is  made  between 
global  and  local  behavior  of  structures.  Only  the  global 
model  is  coupled  to  the  external  flow.  The  loads  gen¬ 
erated  by  the  global  model  are  then  transferred  to  the 
internal  structures  wherever  needed,  locally. 

The  local  models  handle  details,  which  include  FEM 
of  individual  composite  ply-layups,  rigid  body  mecha¬ 
nisms,  redundant  load  paths,  local  buckling,  crack  prop¬ 
agation,  stress  concentration  due  to  manufacturing  im¬ 
perfections,  embedded  cooling,  and  integrated  smart  ma¬ 
terials  and  actuators.  Local  effects  (e.g.  actuator  loads, 
local  buckling,  composite  tailoring)  are  included  via  lower 
order  models.  The  advent  of  HPC  opens  opportunity  to 
couple  the  important  frequencies  of  local  analysis  directly 
to  global  analysis  using  detailed  modeling.  Performed  ju¬ 
diciously,  based  on  a  fundamental  understanding  of  the 
key  phenomena,  it  can  lead  to  enormous  payoffs  in  design 
innovation  and  cost. 

A  comprehensive  treatise  on  fixed-wing  structural 
analysis  and  design,  as  performed  in  the  Boeing  Corn- 


Figure  9:  The  outboard  wing  cross  section  of  a 
Boeing  737-NG  aircraft;  van  Dam  2002 


Figure  10:  A  Boeing  777  FEM  stress  model  with 
details  of  major  internal  components;  Mohaghegh 
2005 

pany  can  be  found  in  Mohaghegh  [43].  The  philosophy 
of  fail  safe  design  has  led  to  the  well-known  practice  of 
discretely  stiffened  wing  and  fuselage  panels  and  mul¬ 
tiply  redundant  two  and  three-piece  primary  bulkheads 
in  the  construction  of  modern  aircraft  structures.  To¬ 
day,  global  analysis  can  be  performed  using  full  aircraft 
FEM  and  RANS.  The  structural  model  involves  large- 
scale  but  usually  linear  FEM  models.  The  global  loads 
are  then  used  for  internal  stress  calculations  using  de¬ 
tailed  FEM  of  the  type  shown  in  Fig.  10.  Most  often, 
this  step  involves  a  static  analysis.  Since  the  1990s,  all 
primary  structures  of  commercial  airplanes  like  the  B777 
and  A340  are  certified  using  such  FEM  analysis.  A  typ¬ 
ical  example  is  the  wing-body  junction  which  includes 
the  cargo  bay  and  landing  gear  —  a  region  of  complex 
redundant  load  paths  and  structural  discontinuities,  see 
Fig.  11.  The  need  for  analysing  such  structures  led  to  the 
original  development  of  sub-structuring  methods.  Com¬ 
posites  materials  are  used  on  selected  components  like 
the  radomes,  fairings,  all  control  surface  panels,  engine 
nacelles,  and  torque  boxes.  The  main  purpose  is  to  pro¬ 
vide  higher  impact,  fatigue,  and  corrosion  resistance. 

Almost  every  modern  rotor  blade  is  an  all  compos¬ 
ite  construction  (although  the  main  spar  is  frequently 
of  Titanium  construction).  The  internal  layout  of  rotor 
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Figure  11:  A  routine  Boeing  747  wing-body  interaction  FEM  stress  analysis  model  with  typical  sub¬ 
structures;  Mohaghegh  2005 


blades,  however,  are  much  simpler  compared  to  the  me¬ 
chanical  complexity  of  fixed  wing  sections.  The  control 
mechanism  for  conventional  main  rotors  are  all  essen¬ 
tially  concentrated  at  the  root  end.  In  the  case  of  rotors, 
the  state-of-the-art  in  global  and  local  models  are  ID 
beam  FEM  and  3D  stress  analysis  respectively.  The  ID 
beam  properties  are  first  obtained  by  a  separation  of  the 
3D  problem  into  a  ID  problem  and  a  2D  cross-sectional 
analysis  problem  (see  section  on  rotary  wing).  Analogous 
to  its  aerodynamic  lifting-line  counterpart,  the  theory  is 
not  meant  for  end  stresses  or  stress  concentrations  at 
the  attachment  point.  More  importantly,  the  root  end 
flex-beams  and  torque  tube  structures  of  modern  hin¬ 
geless  and  bearingless  rotors  require  that  the  external 
loads  be  transferred  correctly  to  internal  stresses.  This 
is  currently  performed  by  a  detailed  re-calculation  using 
commercial  FEM  codes  like  ABAQUS,  ANSYS,  etc,  in  a 
static  fashion  (using  the  dynamic  forcing  obtained  from 
the  global  model).  It  is  desirable  that  future  HPC  based 
simulations  make  this  approximate  re-calculation  redun¬ 
dant. 

While  the  1-D  beam  FEM  model  of  a  rotor  blade 
is  comparatively  simpler,  the  models  include  geometric 
nonlinearities  rigorously  that  are  critical  for  rotors.  Ro¬ 
tary  wing  aeroelasticity  is  inherently  nonlinear,  the  non- 
linearities  stem  from  the  centrifugal  and  Coriolis  forces 
due  to  blade  rotation.  The  Coriolis  nonlinearity  is  in  turn 
coupled  with  the  rotor  trim  state.  This  level  of  coupling 
between  CSD  and  vehicle  dynamics  is  absent  in  fixed 
wings,  and  this  has  its  implications  on  the  efficiency  of 
methods  that  couple  rotorcraft  CSD  with  CFD. 

Aeroelastic  Phenomena  and  CFD/CSD  Predic¬ 
tions 

Flutter 

Flutter  analysis  in  subsonic  regimes  are  satisfacto¬ 
rily  performed  by  standard  k  and  p  —  k  methods  using 
doublet-lattice  type  linear  aerodynamic  models.  These 
options  are  available  today  as  part  of  FEM  codes  like 
NASTRAN  and  ASTROS.  In  supersonic  regimes,  meth¬ 


ods  related  to  piston  theory  are  used.  A  description  of  the 
current  status  of  these  methods,  along  with  their  aeroe¬ 
lastic  applications,  can  be  found  in  Yurkovich  [44].  High- 
performance  fighter  aircrafts  are,  however,  flutter  critical 
in  the  transonic  regime  where  unsteady  shock  motions 
of  three  different  types  determine  the  nature  of  energy 
transfer  between  the  fluid  and  the  structure.  A  linear 
structure  is  still  enough,  only  the  aerodynamic  nonlin- 
earities  need  to  be  predicted.  Hence  high-fidelity  RANS 
CFD  approaches  are  sought. 

Unlike  helicopter  blades,  where  flap-lag  flutter  of 
hingeless  rotors  is  determined  by  structural  nonlineari- 
ties  coupled  to  the  vehicle  operating  state,  fixed  wing 
bending-torsion  flutter  is  a  linear  phenomena.  The  effect 
of  vehicle  g-level  is  only  via  the  effect  of  trim  angle  of  at¬ 
tack  on  the  shock  motion.  At  a  given  angle  of  attack,  and 
dynamic  pressure,  a  conceptually  straight-forward  fluid- 
structure  transient  response  can  be  used  to  extract  the 
frequency  and  damping.  Because  classical  logarithmic 
decay  methods  are  sensitive  to  noise,  refined  eigensys- 
tem  realization  algorithms  are  used  to  extract  these  pa¬ 
rameters  for  the  particular  modes  of  interest.  Most  often 
(except  for  free-body  control  surface  flutter)  these  are  the 
first  few  antisymmetric  modes  of  the  aircraft.  Depend¬ 
ing  on  the  convergence  and  divergence  of  response,  gen¬ 
erally  5  to  6  transient  calculations  are  enough  to  bracket 
a  flutter  point  (i.e.  a  frequency,  damping  vs.  dynamic 
pressure  point  corresponding  to  zero  damping  at  a  given 
Mach  number). 

One  of  the  early  CFD  based  flutter  calculation  for 
a  complete  aircraft  was  performed  by  Melville  in  2001 
and  2002  using  Euler  [45]  and  RANS  [46]  coupling,  for  a 
F-16  fighter  aircraft.  The  linear  FEM  structural  model 
used  was  the  same  as  used  by  Denegri  earlier  [47]  for  a 
doublet-lattice  based  flutter  analysis.  Ten  symmetric  and 
ten  anti-symmetric  mode  shapes  were  considered.  An  ex¬ 
tra  mode  was  incorporated  for  the  leading  edge  flap.  The 
flap  was  represented  as  a  damped  first-order  system  that 
relaxed  into  a  commanded  position.  The  FEM  model 
is  shown  in  Fig.  12(a).  It  was  built  from  several  sim¬ 
pler  models  representing  each  of  the  main  components 
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(b)  Fluid-structure  interface  ge¬ 
ometry  representation 


(c)  Measured  and  predicted  flutter  boundary 


Figure  12:  Measured  and  predicted  flutter  onset  boundaries  for  a  F-16  fighter  aircraft  using 
RANS/FEM  modal  coupling;  Angle  of  attack  0=1.5°;  Melville  2001,  2002 


(a)  Predicted  and  measured  torsion  mode  damping  (b)  Predicted  and  measured  bending  /  torsion  frequencies 


Figure  13:  Predicted  and  measured  aeroelastic  parameters  for  a  F-16  fighter  aircraft  at  3000  m  altitude 
using  RANS/full  FEM  coupling;  Farhat  et  al.  2003 


-  the  fuselage,  wings,  leading-edge  flap,  flaperon,  tip 
launcher,  and  tails.  The  fluid-structure  interface  repre¬ 
sentation  was  constructed  by  seperate  extrapolation  of 
each  component.  The  interface  is  shown  in  Fig.  12(b). 
The  curvature  was  provided  by  the  underlying  structural 
components,  extrapolation  to  surface  was  done  linearly. 
The  RANS  domain  (Baldwin-Lomax)  contained  20  to  36 
overset  grids  with  1.9M  to  3.5M  points.  The  details  of 
the  inlet,  exhaust,  ventral  fins  and  under-wing  store  were 
not  included.  The  coupling  method  and  the  moving  grid 
mechanism  (with  its  associated  Geometric  Conservation 


Law)  was  based  on  an  earlier  work  by  Morton  et  al.  [48]. 
A  single  CPU  was  used  for  structures  and  grid  deforma¬ 
tion,  the  fluids  solver  was  run  in  parallel. 

A  perturbation  response  damped  the  three  primary 
symmetric  modes  into  steady  deflections.  The  three 
primary  anti-symmetric  modes  were  lowly  damped  and 
showed  zero  mean  oscillations.  The  conventional  proce¬ 
dure  for  extracting  modal  damping  is  to  use  classical  log¬ 
arithmic  decay  methods.  For  a  full  FEM  model,  real-time 
parameter  identification  methods,  as  used  in  flight  tests, 
are  used.  The  procedure  is  recognized  as  a  challenge  in 
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(a)  A  Boeing  twin-engine  transport  wind-tunnel  model 
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(b)  Transonic  flutter  onset  dynamic  pressures  for  a  hump  mode  and  a  tip  mode 


Figure  14:  Predicted  and  measured  flutter  onset  for  a  twin-engine  Boeing  transport  aircraft  model 
using  RANS/FEM  modal  coupling;  Hong  et  al.  2003,  reproduced  from  Bartels  and  Sayma  2007 


general.  Often,  a  doublet-lattice  based  linear  analysis 
is  performed  a  priori  to  identify  the  critical  modes.  In 
this  particular  instance  for  example,  linear  analysis  re¬ 
sults  from  Denegri  [47]  was  used  to  identify  the  second 
mode  as  the  critical  mode.  This  mode  was  then  tracked 
during  the  CFD/CSD  simulation.  Figure  12(c)  compares 
the  predicted  flutter  boundary  with  flight  test  data. 

Far  hat  and  his  co-researchers  [49,  50]  have  per¬ 
formed  RANS/full  FEM  simulations  of  the  same  aircraft. 
The  structural  model  was  that  shown  earlier  in  Fig.  5.  It 
contained  around  0.1 7M  linear  DOFs  comprised  of  bar, 
beam,  solid,  plate,  shell,  and  composite  elements.  The 
fluid  solver  was  run  on  up  to  24  CPUs  on  an  SGI  Ori¬ 
gin  3200  computer  with  parallel  speed  up  of  21.  Of  the 
total  CPU  time,  60%  was  spent  on  fluids,  38%  in  mesh 
deformation,  and  only  2%  on  structures.  The  authors 
note  that  for  nonlinear  (geometric)  FEM,  the  structures 
require  15%  of  time.  The  total  time  for  one  transient 
calculation  on  24  CPUs  was  2.5  hours.  An  Eigensystem 
Realization  algorithm  was  used  to  extract  the  frequencies 


and  damping  of  the  two  lowest  modes,  see  Fig.  13.  The 
algorithm  requires  only  two  cycles  of  history  as  long  as 
the  sampling  rate  is  500-1000  Hz.  A  typical  wing-section 
analysis  using  2-D  CFD  coupled  to  a  2-DOF  spring-mass 
system  was  also  performed  for  comparison.  Recently  in 
2007,  Lieu  and  Far  hat  [51]  have  used  the  same  test  case  to 
systematically  compare  these  results  with  those  obtained 
using  Proper  Orthogonal  Decomposition  based  reduced 
order  aerodynamic  models. 

A  description  of  the  status  of  RANS  based  flutter 
predictions  in  commercial  aircraft  design  can  be  found  in 
the  review  paper  by  Bartels  and  Sayma  [40] .  A  collabora¬ 
tive  Boeing-NASA  effort  recently  examined  flutter  onset 
on  a  half-span  twin-engine  transport  model  including  the 
aerodynamic  interactions  of  the  wing,  nacelle,  and  strut, 
the  effect  of  wind  tunnel  wall,  and  the  effect  of  static  aero. 
The  details  of  the  work  originally  appeared  in  Hong  et 
al  [52],  and  are  reproduced  here  in  Fig.  14  from  refer¬ 
ence  [40].  The  CFD  code  used  was  CFL3D  (version  6), 
which  coupled  RANS  to  structural  modes.  A  second- 
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(a)  Flight  test 


(b)  1/72-scale  model  test 


Figure  15:  Flow  visualization  of  LEX  vortices;  (a)  An  LEX  vortex  on  a  F/A-18  High  Angle-of-attack 
Research  Vehicle  (HARV),  wing-tip  view  at  a=25°,  /3=-1.4°;  (b)  A  1/72-scale  model  plan  view  at 
o=30°;  Lee  2000 


order  predictor-corrector  based  coupling  was  used  in  this 
study  between  the  fluid  and  structure.  The  model  exhib¬ 
ited  a  ‘hump’  mode,  i.e.  a  lightly  damped  mode  which 
is  unstable  only  on  a  limited  range  of  dynamic  pressure. 
The  analysis  predicted  this  mode.  Predictions  for  an¬ 
other  unstable  mode,  the  tip  mode,  showed  inconsistent 
results.  Predictions  showed  greater  error  when  the  ac¬ 
tual  cruise  shape  was  used  compared  to  the  undeformed 
jig  shape. 

Tail  Buffet 

Fixed-wing  tail  buffet  studies  date  back  to  the  1930s 
and  officially  began  with  a  Junkers  F13  ge  monoplane 
crash  in  Meopham  in  England.  The  official  British  inves¬ 
tigation  attributed  the  tragic  accident  to  horizontal  tail 
buffeting.  Subsequently,  early  work  focused  mostly  on 
horizontal  tail  buffet.  Since  the  late  1970s  and  80s  the 
buffet  problem  has  received  renewed  attention  with  the 
introduction  of  the  F-14,  F-15,  F-16  and  F/A-18  fighters 
in  the  US  Air  Force  and  Navy.  This  time,  the  focus  has 
been  on  vertical  tail  buffet. 

Modern  fighters  are  designed  with  maneuver  capa¬ 
bilities  at  high  angle  of  attack  (AOA),  typically  25°  to 
35°,  sometimes  as  high  as  60°.  To  generate  lift  at  these 
high  AOA,  a  Leading  edge  Extension  (LEX)  is  designed 
on  each  wing  to  create  two  vortices  that  trail  downstream 
over  the  aircraft.  At  the  same  time  small  disturbances 
from  the  fore-body  are  also  reinforced  by  the  two  LEX 
vortices.  The  vertical  tails  (also  called  fins)  which  pro¬ 
vide  directional  stability  and  control,  are  placed  to  take 
full  advantage  of  these  two  vortices.  However,  if  the  vor¬ 


tices  fail  to  traverse  the  adverse  pressure  gradient  on  the 
wing,  and  burst,  they  impinge  on  the  tail  instead  and 
produce  large  vibration  and  consequent  fatigue  damage. 
This  phenomena  is  called  tail  buffet.  Figure  15  shows 
flow  visualization  of  vortical  flows  that  produce  tail  buf¬ 
fet. 

Tail  buffet  without  the  interaction  with  LEX  vortex 
have  also  been  reported.  For  example,  on  the  F-15,  the 
main  wing  vortex  system  generated  buffet  loads.  It  is 
only  for  the  F/A-18  aircraft  that  a  large  volume  of  ex¬ 
perimental  and  analytical  studies  exist  for  a  systematic 
understanding  of  the  problem.  Lee  [54]  in  2000  has  com¬ 
prehensively  described  the  problem,  and  reviewed  the 
status  of  fundamental  understanding,  experimental  data, 
and  analytical  predictions  of  the  tail  buffet.  The  problem 
has  acquired  increased  relevance  as  the  three  largest  U.S. 
fighter  aircraft  programs  have  twin  tails,  F/A-18E/F,  F- 
22,  and  the  F-35  (Joint  Strike  Fighter). 

There  have  been  several  noteworthy  attempts  at  al¬ 
leviating  tail  buffet  from  the  late  80s  through  90s  using 
extensive  wind  tunnel  and  flight  tests.  The  only  method 
used  today  is  the  LEX  fence  developed  by  McDonnell 
Douglas  Corporation  via  trial  and  error  using  wind  tun¬ 
nel  experiments  and  flight  tests.  The  LEX  fence  is  a 
flat  trapezoidal  plate  standing  on  the  LEX  upper  surface 
aligned  in  a  stream-wise  direction. 

The  complexity  of  the  problem  and  its  sensitivity 
to  flight  conditions  and  aircraft  configurations  makes  it 
ideally  suited  for  a  high-fidelity  analysis.  At  the  time  of 
Lee’s  review  there  was  only  limited  RANS  calculation  of 
the  tail  buffet  problem.  These  early  calculations  were  by 
Rizk  and  Gee  [55]  in  1992,  Ghaffari  et  al  [56]  in  1993,  and 
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(a)  Instantaneous  streamlines  of  LEX  and  wing 
vortices  without  LEX  fence 


(b)  Instantaneous  streamlines  of  LEX  and  wing 
vortices  with  LEX  fence 
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(c)  Without  LEX  fence 


(d)  With  LEX  fence 


Figure  16:  Predicted  and  measured  tail  bending  moments  during  F/A-18  tail  buffet  alleviation  us¬ 
ing  LEX  fences;  RANS/FEM  coupling,  FEM  of  vertical  tail  only;  Bending  moments  are  integrated 
differential  pressure  moments,  not  structural  gauge  loads;  Sheta  2004 


Gee  et  al.  [57]  in  1996.  These  studies  focussed  on  CFD, 
the  structure  was  assumed  to  be  rigid,  and  had  limited 
success  in  predicting  structural  buffet  loads  on  the  tail. 

A  RANS/FEM  prediction  of  buffet  loads  on  the 
F/A-18  tail  was  carried  out  by  Sheta  in  2004  [58].  The 
study  clarified  the  fundamental  contribution  of  the  LEX 
fence  in  alleviating  buffet  loads.  The  aircraft  was  as¬ 
sumed  to  be  rigid,  but  the  tail  was  modeled  using  3- 
D  hexahedral  brick  FEM  (isotropic,  Aluminum).  The 
RANS  domain  had  53  structured  blocks  with  a  total  of 
2.68  M  grid  points.  The  total  time  required  on  6  CPUs 
(of  1.0  GHz  speed  each)  was  142  h.  Figure  16  shows  the 
predicted  flowfields  and  tail  bending  moments  from  this 
study  with  and  without  the  LEX  fences.  The  off-body 
flow  was  considered  to  be  laminar  in  this  study. 

Recently,  Morton  et  al.  [59]  have  used  detached- 


eddy  simulation  (DES)  based  turbulence  model  for  a 
more  meaningful  calculation  of  vortex  breakdown.  DES 
is  a  popular  RANS/LES  hybrid  model.  It  uses  RANS 
in  the  boundary  layer  and  LES  in  the  regions  of  mas¬ 
sive  separation.  Designed  primarily  for  external  flows, 
it  can  be  implemented  by  a  relatively  simple  modifica¬ 
tion  of  an  existing  Spalart-Allmaras  (S-A)  RANS  model 
in  an  unsteady  solver.  The  LES  domain  does  not  re¬ 
quire  an  explicitly  defined  SGS  model.  The  modification 
to  the  S-A  model  automatically  implies  the  existence  of 
a  Smagorinsky-type  SGS  model  under  the  assumption  of 
local  equilibrium  between  the  production  and  dissipation 
of  eddy  viscosity.  The  simulations  were  performed  for  the 
F-18  High  Angle-of- Attack  Research  Vehicle  (HARV). 
For  the  chosen  flight  condition,  the  F-18  HARV  configu¬ 
ration  matches  an  F/A-18C  configuration  with  leading- 
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(a)  LEX  vortex  breakdown  around  a 
F/A-18  aircraft  using  DES  model 
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(b)  Vortex  breakdown  location  (nose  is  x=0) 
validated  with  flight  test  data 
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Figure  17:  RANS  calculation  of  vortex  bursting 
and  tail  impingement  using  Detached  Eddy  Sim¬ 
ulation  (DES);  no  LEX  fences;  Morton  et  al.  2007 


edge  flaps  set  to  -33°,  trailing-edge  flaps  set  to  0°,  and 
the  diverter  slot  through  the  LEX  fence  closed.  The  an¬ 
gle  of  attack  was  30°,  Mach  number  0.2755,  Reynolds 
number  13M  (corresponds  to  altitude  of  20,000  ft).  The 
DES  simulations  were  performed  on  a  baseline  unstruc¬ 
tured  grid  of  3.6M  cells  (using  the  commercial  CFD  code 
Cobalt).  A  time-averaged  solution  was  used  to  produce 
an  AMR  grid  with  3.9M  cells  using  Pirzadeh’s  grid  adap¬ 
tation  method  [60].  The  analysis,  however,  did  not  in¬ 
clude  a  structural  dynamic  model;  the  focus  was  on  the 
flow  field. 

Limit  Cycle  Oscillations  (LCO) 

Fundamental  understanding  of  LCO  and  its  satis¬ 
factory  prediction  remain  beyond  the  state  of  the  art  in 
fixed-wing  aeroelasticity.  The  LCO  is  a  general  term.  In 


fixed  wing  fighters  it  is  characterized  by  sustained,  pe¬ 
riodic,  antisymmetric  oscillations  of  the  wing  which  are 
not  catastrophically  divergent.  Inside  the  fuselage,  it  is 
felt  as  a  lateral  motion.  It  has  been  a  persistent  prob¬ 
lem  since  the  mid-1970s  in  high  subsonic  and  transonic 
speeds  for  configurations  with  missiles  on  the  wingtips 
and  heavy  stores  on  the  outboard  pylons.  The  F-16  and 
F/A-18  fighters  with  wing  tip  missile  launchers  tradition¬ 
ally  encountered  LCO.  Similar  phenomena  are  expected 
for  the  F-22  and  F-35  as  flight  tests  begin. 

Bunton  and  Denegri  [61]  presents  an  insightful  com¬ 
mentary  on  fixed-wing  LCO.  The  phenomenon  is  similar 
to  classical  flutter,  in  that  the  oscillations  occur  at  a  sin¬ 
gle  frequency  close  to  that  predicted  by  linear  flutter. 
The  mode  shape,  at  the  onset,  also  resembles  that  pre¬ 
dicted  by  linear  flutter.  In  level  flight,  the  onset  speed 
lie  quite  close  to  that  predicted  by  linear  flutter.  As  a  re¬ 
sult  LCO  has  been  termed  limit  cycle  flutter  or  constant 
amplitude  flutter  by  many,  a  terminology  not  preferred 
by  others  as  the  word  ‘flutter’  has  historically  referred  to 
catastrophic  divergence.  Many  speculate  that  the  gene¬ 
sis  of  LCO  is  not  the  classical  bending-torsion  flutter  at 
all,  but  the  aircraft  flight  control  system. 

During  LCO,  the  amplitude  remains  constant  for  a 
given,  stable  flight  condition.  If  the  aircraft  is  accelerated 
to  a  higher  speed,  the  amplitude  would  increase  and  set¬ 
tle  into  a  higher  level  once  the  higher  speed  is  reached.  If 
an  attempt  is  made  to  back  out  of  LCO,  the  oscillations 
often  persist  below  the  original  onset  values  (hysteresis). 
In  addition  to  LCOs  involving  the  entire  aircraft-wing 
store  structure,  LCOs  can  occur  involving  localized  non- 
linearities  like  control  surface  free-play. 

LCO  occurs  both  in  level  flight  as  well  as  in  high 
load  factor  maneuvers.  Usually  the  onset  speed  is  lower 
in  the  latter.  Compared  to  LCO  onset  in  level  flight, 
those  occurring  in  maneuvers  are  much  less  predictable. 
For  example,  LCO  magnitudes  can  increase  from  lg  to  4g 
and  then  begin  decreasing  so  as  to  disappear  at  5.5g  [61]. 
For  a  different  store  configuration  there  may  not  be  any 
LCO  at  all  until  5g,  at  which  LCO  begins,  and  grow 
so  rapidly  so  as  to  resemble  the  onset  of  classical  catas¬ 
trophic  flutter.  The  key  unresolved  challenge,  beyond 
that  of  predicting  the  onset,  is  predicting  the  magnitude 
of  LCO.  In  fixed  wing,  unlike  rotorcraft,  the  fatigue  of 
components  is  of  lesser  concern  compared  to  the  limit 
loads.  (Fatigue  concerns  are  evident  however  for  test  air¬ 
craft  repeatedly  exposed  to  flutter  tests  and  LCO.)  LCO 
concerns  are  related  to  weapons  accuracy  for  smart  and 
unguided  weapons,  and  reliability  of  ordnance  release. 
The  transonic  aerodynamic  nonlinearity  that  can  lead  to 
LCO  has  been  explained  by  Bendiksen  [62]  based  on  the 
categorization  of  shock  motions  by  Tij  deman  and  Seebass 
in  1980  [63].  It  was  shown  that,  for  a  clean  wing,  when 
a  bending-torsion  flutter  mode  triggers  a  transition  from 
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(a)  External  stores  on  a  F-16  fighter  during  flutter  and  LCO 
flight  tests 


Mach  Number,  Mt 


(b)  Measured  maximum  wing-tip  LCO  response  in 
level  flight 

Figure  18:  Measured  wing-tip  LCO  deformations 
of  a  F-16  fighter  aircraft  with  under-wing  stores; 
Denegri  et  al.  2005 


(a)  Configuration  1 


Mach  Number,  Mm 

(b)  Configuration  2 

Figure  19:  Predicted  and  measured  LCO  magni¬ 
tudes  using  RANS  CFD/FEM  coupling  in  level 
flight;  (a)  Effect  of  slendor  body- wing  model  for 
launchers,  (b)  LCO  predicted  only  with  launcher 
model;  Thomas  et  al.  2007 


type-A  (shock  motion  oscillates  sinusoidally)  to  type-B 
(shock  disappears  during  a  downstream  motion)  shock 
motion,  then  the  associated  rapid  decrease  of  aerody¬ 
namic  work  done  per  cycle  can  lead  to  LCO.  The  role 
played  by  under-wing  stores  and  structural  nonlineari- 
ties  including  free-play  are  more  configuration  specific. 
As  flight  testing  becomes  more  and  more  expensive,  and 
the  number  and  combinations  of  under-wing  store  con¬ 
figurations  increase,  higher-fidelity  methods  are  ideally 
suited  for  a  comprehensive  attack  on  the  LCO  problem. 

One  of  the  early  high-fidelity  attempt  to  simulate 
LCO  was  that  of  Melville  in  2002  [46]  described  earlier 
in  the  section  on  flutter.  While  calculating  the  onset  of 
flutter  he  investigated  whether  the  undamped  response 
(or  lowly  damped  response  near  the  flutter  boundary) 
resembled  an  LCO  as  seen  in  flight.  What  distinguishes 
an  LCO  response  from  a  zero-damping  periodic  flutter 


point  is  that  the  LCO  amplitude  is  independent  of  the 
initial  excitation,  and  depends  only  on  the  underlying 
limiting  mechanism.  The  second  mode  excitation,  when 
calculated  with  different  set  of  initial  excitation,  settled 
to  a  different  amplitude  (proportional  to  the  initial  exci¬ 
tation).  Thus  it  was  not  LCO. 

More  recently,  in  2005,  Denegri  et  al.  [64]  have  gener¬ 
ated  benchmark  LCO  data  for  the  F-16  aircraft.  A  sam¬ 
ple  of  benchmark  measured  LCO  deflections  along  with 
the  associated  wing  measurement  locations  are  shown  in 
Fig.  18.  Oscillation  frequencies  range  from  7.9Hz  to  8.2 
Hz  with  lower  frequencies  around  Mach  0.9. 

In  2007,  Thomas  et  al.  [65]  have  used  a  nonlinear 
harmonic  balance  based  RANS  solver  to  calculate  the 
LCO  magnitudes  for  different  configurations  of  the  F-16. 
A  slender-body  wing-theory  was  used  to  model  the  ex- 


17 


ternal  weapons  and  wing  stores  and  study  their  impact 
on  LCO  prediction  for  eight  different  weapons  and  store 
configurations.  A  comparison  of  the  predicted  LCO  mag¬ 
nitudes  with  those  of  measured  data,  as  in  Fig.  18(b),  is 
given  in  Fig.  19(a).  It  corresponds  to  2000  ft  altitude, 
1.5°  angle  of  attack  flight  test.  The  figure  is  sufficient 
to  highlight  the  poor  state  of  affairs  in  LCO  prediction. 
Figure  19(b)  compares  predictions  for  a  different  config¬ 
uration,  but  at  the  same  altitude  and  angle  of  attack. 
The  LCO  could  at  least  be  predicted  with  the  wing  store 
model. 

Parker  et  al.  [66]  have  carried  out  systematic  studies 
on  the  effect  of  under-wing  stores  on  LCO  using  a  simple 
rectangular  cantilevered  wing  model.  The  details  of  the 
wing  model  can  be  found  in  [67].  The  model  is  based 
on  Eastep  and  Olsen’s  heavy  version  [68]  of  the  origi¬ 
nal  Goland  wing  [69].  They  found  that  the  role  of  store 
aerodynamics  was  to  transfer  additional  energy  into  the 
structure  increasing  the  magnitude  of  limit  cycle  oscilla¬ 
tion. 

In  addition  to  the  F-16  and  F/A-18,  LCO  studies 
have  been  documented  on  a  B-l  like  configuration  [71] 
and  on  a  B-2  bomber.  The  latter  phenomena  is  now  well- 
known  as  the  Residual  Pitch  Oscillation  (RPO)  and  has 
been  studied  by  Dreim  and  his  co-researchers  [72]  using 
the  Computational  Aeroelasticity  Program  -  Transonic 
Small  Disturbance  (CAP-TSDv)  code. 


Figure  20:  FEM  model  of  bladed  disks:  (a)  a  disk 
sector;  (b)  blade  contact  surface;  (c)  disk  contact 
surface;  and  (d)  an  area  contact  element;  Petrov 
and  Ewins  2006. 


TURBOMACHINES 

The  advances  in  fundamental  understanding,  analy¬ 
sis,  and  prediction  of  turbomachinery  aeroelasticity  can 
be  broadly  traced  from  1970s  to  the  present  time  by  sur¬ 
vey  papers  in  Mikolajczak  et  al.  in  1975  [73],  Bendikson 
in  1990  [74],  Marshall  and  Imregun  in  1996  [75],  and  Bar¬ 
tels  and  Sayma  in  2007  [40].  The  historical  evolution  of 
the  aerodynamic  methods  has  been  traced  recently  by 
Cumpsty  and  Greitzer  [76].  Structural  dynamics  (and 
flutter)  has  been  surveyed  by  Srinivasan  in  1997  [77], 
Slater  in  1999  [78],  and  more  recently  Castanier  and 
Pierre  in  2006  [79]. 

The  purpose  of  the  following  sections  is  to  review 
the  status  of  turbomachinery  high-fidelity  computational 
aeroelasticity.  As  in  the  case  of  fixed  wing,  the  focus  here 
is  again  on  CSD,  and  coupled  high-fidelity  CFD/CSD 
methods.  The  first  section  summarizes  the  status  of 
detailed  CSD.  The  second  section  summarizes  the  key 
aeroelastic  phenomena,  and  the  status  of  CFD/CSD  in 
predicting  them. 

Structural  Modeling 

Finite  element  models  of  turbine  blades  comprise  of 
millions  of  DOFs  made  up  to  3-D  linear  solid  elements. 
Commercial  FEM  packages  like  ANSYS  and  ABAQUS 
are  routinely  used  for  industrial  as  well  as  research  pur¬ 
poses.  Turbine  blades  are  made  out  of  nickel  based  single¬ 
crystal  super  alloys.  In  addition  to  unsteady  stresses, 
turbine  blades  have  to  endure  high  temperature  (par¬ 
ticularly  during  start-up  and  shutdown)  and  exposure 
to  high-pressure  hydrogen  (in  case  of  internal  coolant). 
Furthermore,  the  blades  are  prone  to  severe  fretting  by 
the  root  end  friction  damper,  and  contact  damage  at  the 
dovetail  attachments  which  can  lead  to  crystallographic 
initiation  and  crack  growth.  Detailed  3D  FEM  analyses 
are  carried  out  with  special  surface-to-surface  contact  el¬ 
ements  for  such  purposes,  see  for  example  Arakere  et 
al.  [80]  and  the  references  therein. 

The  forced  response  of  bladed  disks  are  affected  by 
contact  surfaces  at  the  blade-disk  dovetail  attachments, 
part-span  interlocking  shrouds,  and  under-platform  fric¬ 
tion  dampers  at  the  root.  These  components  are  used  to 
avoid  dangerous  resonance  regimes.  Significant  amount 
of  research  is  conducted  to  model  these  discontinuity  phe¬ 
nomena  within  3-D  FEM.  Recent  work  on  3-D  FEM  mod¬ 
eling  of  dovetail  attachments  can  be  found  in  Beisheim 
and  Sinclair  [81].  Recent  work  on  modeling  of  under¬ 
platform  friction  dampers  can  be  found  in  Petrov  et 
al.  [82,  83].  Figure  20,  taken  from  the  latter  work,  illus¬ 
trates  the  placement  of  an  area  contact  element  within  a 
typical  FEM  model  of  a  bladed  disk. 

When  all  sectors  of  a  bladed  disk  are  similar,  the 
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(a)  Three  nodal  diameter  mode  of  a  tuned  bladed  disk 

Figure  21:  Mode  shapes  of  a  tuned  and  mistuned 

structure  is  called  tuned,  or  having  cyclic  symmetry.  The 
analysis  can  then  be  simplified  to  one  or  a  few  sections. 
Mistuning  destroys  cyclic  symmetry  and  produces  mode 
localization  in  a  bladed  disk.  This  phenomena  is  partic¬ 
ularly  important  for  High  Cycle  Fatigue  (HCF),  where 
high  localized  stresses  can  lead  to  failure.  Figure  21 
shows  an  example  of  a  modern  bladed  disk  FEM  model 
showing  the  effect  of  mistuning  on  mode  shapes.  Studies 
on  mistuning  have  received  attention  since  1970s,  led  by 
the  pioneering  work  of  Whitehead  [84]  and  Ewins  [85]. 
However,  it  is  only  recently  that  the  emergence  of  whole 
annulus  multi-row  CFD  have  made  reliable  forced  re¬ 
sponse  predictions  of  mistuned  bladed  disks  appear  fea¬ 
sible  in  near  future.  For  a  recent  review  of  the  vari¬ 
ous  modeling  and  analysis  techniques,  and  the  assess¬ 
ment  of  mistuning  sensitivity,  see  Ref.  [79].  The  mode 
localization  due  to  mistuning,  as  mentioned  earlier,  is 
more  pronounced  in  the  presence  of  frequency  veering. 
Frequency  veering  is  a  common  phenomenon  in  rotat¬ 
ing  structures.  Coupled  modes  approach  each  other  with 
increase  in  RPM  and  veer  away  after  exchanging  mode 
shapes  due  to  changing  contributions  of  centrifugal  load¬ 
ing  in  each.  The  same  phenomena  occurs  in  helicopter 
blades,  except  that  the  low  aspect  ratio  of  the  turbine 
blades  (around  1.5-3  compared  to  10-15  for  rotorcraft) 
and  the  large  flexible  disk  compress  a  larger  number  of 
vibration  modes  within  a  smaller  frequency  range  which 
are  richer  in  interaction.  A  unique  phenomena,  not  to 
be  found  in  rotorcraft  blades,  is  the  interaction  between 
disk  dominated  and  blade  dominated  modes.  It  is  ob- 


(b)  Localized  mode  shape  of  a  mistuned  bladed  disk 


modern  bladed  disk;  Castanier  and  Pierre  2006 

served,  that  the  sensitivity  of  a  bladed  disk  to  mistuning 
increases  with  frequency  veering.  Recent  work  on  mis¬ 
tuning  and  mode  localization  can  be  found  in  Petrov  and 
Ewins  [86],  Rivas-Guerra  and  Mignolet  [87],  and  Kenyon 
et  al.  [88,  89]. 

The  whole  annulus  treatment  of  a  bladed  disk,  with 
or  without  mistuning,  is  in  itself  an  idealization  which 
holds  true  only  for  isolated,  dynamically  independent 
rotors.  A  disk  with  identical  blades  is  not  necessar¬ 
ily  tuned,  particularly  for  multi-stage  compressor  blades. 
The  modes  of  one  stage  interact  with  those  of  its  neigh¬ 
bor,  which  usually  contains  a  different  number  of  blades 
giving  rise  to  mistuning  of  the  coupled  stages.  Unlike 
helicopters  (co-axial)  the  stage  to  stage  spacing  is  only  a 
fraction  of  the  blade  chord.  In  addition,  the  disk  is  large 
and  flexible,  compared  to  the  hub  in  helicopters.  Bladh 
et  al.  [90]  have  studied  the  effects  of  structural  coupling 
of  multi-staged  bladed  disks,  see  Fig.  22. 

Aeroelastic  Phenomena  and  CFD/CSD  Predic¬ 
tions 

Aeroelastic  Phenomena 

A  modern  turbofan  engine  and  a  schematic  illustrat¬ 
ing  its  instability  prone  components  are  shown  in  Fig.  23. 

Static  aeroelasticity  in  turbomachines  deals  with  the 
calculation  of  blade  un-twist  and  un-camber  because  of 
high  centrifugal  and  pressure  loads.  There  is  no  diver¬ 
gence  phenomenon.  Typical  un-twists  can  be  around  5° 
and  changes  the  stagger  angle.  The  manufactured  shape 
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(a)  FEM  for  single  and  two-stage  rotor  models  (b)  stage  1  forced  response  from  Engine  Order  10  excitation, 

using  tuned  and  mistuned  FEM  single  and  two-stage  models 


Figure  22:  Effect  of  multistage  rotor  structural  coupling  in  turbomachines;  Bladh  et  al.  2003 
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(a)  The  GE  90-115  B  growth  engine  of  the  Boeing  777 


(b)  Schematic  of  a  modern  turbofan 


Figure  23:  The  components  of  a  modern,  high  bypass  ratio  turbofan  engine 


(cold)  is  calculated  based  on  the  anticipated  un-twist  and 
un-camber  of  the  operating  shape  (running).  The  effect 
is  pre-dominant  in  fan  blades.  A  recent  study  can  be 
found  in  Wilson  et  al.  [91]. 

Forced  response  is  most  critical  for  turbine  blades. 
The  High  Cycle  Fatigue  phenomenon  is  a  major  issue 
is  gas  turbine  blades.  In  addition  to  unsteady  forcing, 
high  thermal  stresses  and  internal  coolant  pressure  lead 
to  fatigue  failure.  The  U.S.  Air  Force  estimates  that  55% 
of  jet  engine  safety  Class  A  mishaps  (over  $1  million  in 


damage  or  loss  of  aircraft)  and  30%  of  all  jet  engine  main¬ 
tenance  costs  are  due  to  HCF  [92].  The  forced  response 
phenomenon  also  occurs  in  compressor  blades.  But  here, 
it  occurs  not  only  at  Engine  Order  (EO)  excitations  (sim¬ 
ilar  to  blade  passage  frequency  in  helicopter  main  rotors) , 
but  also  at  Low  Engine  Order  excitations  (LEO).  The  lat¬ 
ter  arise  due  to  stage-to- stage  interactions.  Compared  to 
turbines,  compressors  typically  have  greater  number  of 
blade  stages  which  produce  a  greater  number  of  engine 
order  excitation  frequencies. 
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Figure  24:  A  typical  compressor  map  showing  the 
various  flutter  boundaries 


Flutter  occurs  in  fan  blades,  the  front  stages  of  axial 
compressor  blades,  and  occasionally  in  the  low  pressure 
turbine  blades.  Fan  blades  are  the  most  flexible  compo¬ 
nents,  especially  those  in  the  high  bypass  ratio  commer¬ 
cial  airplane  turbofan  engines.  In  case  of  compressors, 
demands  of  increasing  pressure  ratios  and  weight  reduc¬ 
tion  result  in  highly  loaded  stages  operating  in  transonic 
inflow.  In  addition,  the  front  stages  are  subjected  to  se¬ 
vere  inlet  flow  distortions. 

The  nature  of  flutter  is  complex  in  turbomachines, 
and  the  fundamental  mechanisms  are  often  difficult  to 
distinguish.  Conventionally,  the  various  flutter  regions 
of  a  compressor  are  demarcated  as  shown  in  Fig.  24. 
The  mass  flow  is  a  corrected  value,  i.e.  an  equivalent 
mass  flow  under  ambient  conditions.  The  pressure  ratio 
is  the  ratio  of  exit  stagnation  pressure  to  inlet  stagna¬ 
tion  pressure  of  the  compressor.  The  surge  line  is  the 
highest  operating  line  above  which  the  airflow  can  re¬ 
verse  abruptly  creating  rotating  stall,  called  surge.  The 
subsonic  and  supersonic  stall  flutter  are  two  of  the  more 
distinct  regimes.  The  subsonic  stall  flutter  occurs  in  part- 
speed  conditions  when  the  inlet  flow  is  high  subsonic  or 
transonic.  Typically,  they  are  a  Single  Degree  of  Free¬ 
dom  (SDOF)  torsion  mode  flutter.  The  supersonic  stall 
flutter  occurs  at  high-speed  operating  conditions  when 
the  inlet  flow  is  supersonic  and  contains  detached  lead¬ 
ing  edge  shocks.  Often,  they  are  a  SDOF  bending  mode 
flutter.  A  SDOF  flutter,  unlike  a  classical  frequency  coa¬ 
lescence  type  bending-torsion  flutter  involves  energy  ex¬ 
change  between  a  single  mode  and  the  fluid.  The  driving 
phenomena  in  compressors  are  either  oscillating  shock 
waves  or  stalled  flow.  A  simple  linear  structural  model  is 
enough  as  long  as  the  aerodynamic  damping  is  correctly 
evaluated.  The  complexity  mainly  lies  in  the  latter  task. 


CFD/CSD  Predictions 

The  1990s  saw  the  development  of  large  scale  paral¬ 
lel  solvers  for  unsteady  aerodynamics  in  turbomachines, 
see  for  example  Refs.  [93,  94,  95].  The  high  computa¬ 
tional  requirements  restricted  the  simulations  to  a  single 
stage.  More  recently,  solvers  with  massive  parallelism, 
i.e.  with  demonstrated  scalability  on  100s  of  CPUs  have 
been  developed  [96,  97].  For  example,  the  Stanford  Uni¬ 
versity  code  TFL02000  has  been  tested  for  scalability  on 
1024  processors  for  a  model  problem  with  9  whole  annu¬ 
lus  blade  rows  with  93. 8M  grid  points.  There  have  been 
limited  efforts  however,  to  couple  it  subsequently  to  tur¬ 
bomachinery  CSD,  see  Doi  [98].  Successful  coupling  with 
CSD  was  carried  out  later  at  the  University  of  Maryland 
for  helicopter  main  rotor  loads  and  acoustics  prediction 
(see  section  on  rotorcraft  aeroelasticity) . 

Researchers  have  tried  applying  large  scale  NS  CFD, 
i.e.  whole-annulus,  multi-blade  row,  unsteady  simula¬ 
tions,  to  forced  response  and  flutter  since  the  early  2000s. 
These  were  uncoupled  simulations  —  either  predefined 
blade  motions  were  prescribed  to  CFD  to  study  the  flow 
distortions  [99,  100,  101],  or  CFD  calculated  flow  distor¬ 
tions  were  imposed  on  the  structure  to  calculate  blade 
deformations  [102,  103].  The  first  procedure  can  be  used 
for  an  approximate  determination  of  flutter.  The  second 
procedure  can  be  used  for  an  approximate  determina¬ 
tion  of  forced  response.  The  first  procedure  exploits  the 
assumption  that  flutter  in  turbomachinery  is  predomi¬ 
nantly  a  SDOF  phenomena.  Thus  prescribing  the  torsion 
mode  (or  bending  mode  for  supersonic  flutter)  at  its  nat¬ 
ural  frequency,  and  then  integrating  blade  airloads  over 
each  cycle  to  calculate  damping  is  accepted  as  a  reason¬ 
able  approach. 

Pioneering  work  on  coupled  analysis  began  during 
the  1990s  used  frequency  domain  coupling,  see  for  exam¬ 
ple  Williams  et  al.  [104]  for  panel  solvers,  and  Geroly- 
mos  [105]  for  Euler  solvers.  In  the  work  by  Geroly- 
mos,  the  structural  dynamics  (for  a  particular  mode)  was 
solved  in  the  frequency  domain.  The  periodic  response 
was  then  transferred  to  the  3-D  Euler  solver.  The  Eu¬ 
ler  solver  was  marched  in  time  domain,  until  periodicity. 
This  cycle  of  periodic  CSD  and  time  domain  CSD  was 
repeated  until  convergence.  Note  that  this  approach,  if 
applied  to  a  helicopter  rotor,  will  lead  to  rapid  diver¬ 
gence.  A  similar  method,  but  with  a  key  innovation,  is 
used  in  helicopter  rotors.  This  method,  termed  the  delta 
method,  is  discussed  later  in  the  section  on  rotary  wing 
aeroelasticity. 

Time  domain  coupling  using  staggered  algorithms 
began  to  be  applied  in  the  2000s  with  single  sector  cal¬ 
culations  for  a  single  blade-row.  Carsten  et  al.  [106]  have 
systematically  compared  predicted  modal  decay  between 
the  time  domain  and  frequency  domain  solvers  in  the 
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(a)  The  first  bending  and  first  torsion  eigen- 

modes  (b)  Comparison  of  logarithmic  decrement  for  frequency  and  time  domain 

method,  first  bending(left)  and  first  torsion(right) 


Figure  25:  RANS/FEM  flutter  calculations  for  a  4-bladed  segment  of  a  low  pressure  compressor  annulus 
using  5  structural  modes  (20  modes  in  total);  8  Pentium-4  CPUs;  Carsten  et  al.  2003 


(a)  Comparison  of  two  fan  intake  geometries 


(b)  Least  stable  mode  vs  %  speed 


Figure  26:  A  RANS/FEM  study  of  intake  duct  effects  on  fan  flutter  stability;  Vahdati  et  al.  2002 


presence  of  oscillating  shocks.  For  the  conditions  studied 
no  significant  differences  were  found,  as  shown  in  Fig.  25. 
Whole  annulus  calculations  have  been  performed  exten¬ 
sively  by  researchers  at  the  Imperial  College,  London. 
The  original  work  by  Imregun  et  al.  [107,  108]  showed 
that  forced  vibration  predictions  remained  unchanged 
between  a  typical  sector  analysis  (even  with  adjustment 
of  blade  numbers)  and  whole  annulus  calculations,  as 
long  as  multi  blade  row  interactions  were  not  impor¬ 
tant.  Subsequent  work  by  Sayma  et  al.  [109],  Breard 
et  al.  [110],  and  Vahdati  et  al.  [Ill]  focused  on  multi 
blade  row  simulations.  The  last  study  showed  the  im¬ 
pact  of  intake  duct  geometry  on  the  flutter  stability  of 
inlet  fans,  see  Fig.  26.  Flutter  and  forced  response  calcu¬ 
lations  for  core  compressors  are  more  involved,  compared 
to  fan  flutter,  because  of  the  close  spacing  of  blade  rows. 


Rotor-stator  blade  row  interactions  involve  multiple  up¬ 
stream  and  downstream  rows.  In  2005,  Wu  et  al  [112] 
demonstrated  the  initial  feasibility  of  a  massive  RANS 
simulation  (Baldwin-Barth  turbulence  model)  including 
17  blade  rows,  see  Fig.  27(a).  It  was  performed  on  128 
CPUs  on  SGI  Origin  3000  with  a  relatively  coarse  (com¬ 
pared  to  the  expected  requirements)  grid  of  68M  points. 
A  more  detailed  study  on  a  5  blade  row  model  was  car¬ 
ried  out  recently  by  Vahdati  et  al.  [113],  see  Figs.  27(b) 
and  27(c).  The  structural  model  in  all  of  the  above  stud¬ 
ies  was  elementary  in  comparison  to  CFD.  Only  rotors 
1  and  2  were  modeled,  and  only  the  1st  flap  and  1st 
torsion  modes  (of  various  nodal  diameter  families)  were 
considered. 

Fleeter  and  his  co-researchers  have  performed  fully 
coupled  FEM  based  CFD/CSD  simulations  of  turboma- 
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(b)  (c) 

Figure  27:  A  HPC  based  Whole-annulus,  multi  blade  row  aeroelastic  analysis  of  axial-flow  compressors 
using  unstructured  RANS  /  FEM  coupling;  (a)  A  core-compressor  coarse  17  blade  row  mesh;  (b)  A 
refined  5  blade  row  mesh;  (c)  Structural  mode  11,  susceptible  to  32nd  engine  order  excitation;  Wu  et 
al.  2005;  Vahdeti  et  al.  2007 


chinery  flows  both  for  flutter  [114]  and  for  forced  re¬ 
sponse  [115].  The  latter  work  involved  the  vibratory 
stress  predictions  of  a  1.5-stage  axial  flow  compressor 
(Purdue  Transonic  Compressor)  for  a  inlet  guide  vanes- 
rotor-stator  configuration  of  18-19-18  blades.  Only  the 
rotor  and  stator  were  modeled.  One  segment  from  each 
blade  row  was  placed  in  one  CPU.  The  steel  stator  vanes 
were  modeled  using  2464  hexahedral  brick  finite  ele¬ 
ments,  4  across  chord  and  14  across  span.  The  third 
mode  frequency  was  5203  Hz,  close  to  the  part-speed 
blade  passage  frequency  of  4750  Hz  (part-speed  operat¬ 
ing  RPM  is  15,000,  number  of  blades  19).  A  simulation 
was  performed  at  5203  blade  passage  frequency  (with 
RPM  16,430).  First,  the  unsteady  periodic  flow  of  the 
rotor  was  calculated.  This  required  31  rotor  revolutions 
(80h  of  CPU  time  on  2  IBM  SP2  processors).  The  stator 
vanes  were  not  allowed  to  deform  during  this  time.  This 
was  an  aerodynamic  only  simulation  for  which  a  large 
time  step  (3844  per  period)  could  be  taken.  Next,  the 
fully  coupled  transient  fluid-structure  simulation  was  ini¬ 
tiated.  For  this  case,  the  time  step  was  reduced  by  three 
times.  The  solution  was  marched  until  steady  state  pe¬ 
riodic  conditions.  This  required  700h  of  CPU  time  on  2 
IBM-SP2  processors. 

ROTARY  WING  AIRCRAFT 

A  recent  review  of  rotary  wing  CFD/CSD,  detailing 
the  historical  evolution  and  current  status  of  airloads  and 


structural  loads  prediction  as  of  2006,  can  be  found  in 
Datta  et  al.  [116].  The  following  sections  present  a  sum¬ 
mary  of  that  review.  In  addition,  the  reader  is  brought 
up  to  date  with  the  latest  developments. 

As  before,  the  focus  here  is  on  CSD,  and  coupled 
high-fidelity  CFD/CSD  methods.  A  thorough  review  of 
rotorcraft  CFD  methods,  from  its  early  inception  to  its 
evolution  into  its  present  status  to  future  challenges,  can 
be  found  in  Strawn  et  al.  [117]  in  2007.  As  before,  the 
review  is  divided  into  two  sections.  The  first  section  sum¬ 
marizes  the  status  of  detailed  CSD.  The  second  section 
summarizes  the  key  aeroelastic  phenomena,  and  the  sta¬ 
tus  of  CFD/CSD  in  predicting  them. 

Structural  Modeling 

The  non-linear  coupled  PDEs  of  rotating  blade 
dynamics  were  systematically  derived  in  the  1970s 
by  several  researchers,  notably  Hodges  and  his  co¬ 
researchers  [118,  119],  Kvaternik  et  al.  [120],  and  Rosen 
and  Friedmann  [121].  The  rotor  blade  was  treated  as 
isotropic  and  straight  Euler-Bernoulli  beams.  By  retain¬ 
ing  the  second-order  effects  of  elastic  motion  in  the  ki¬ 
netic  and  strain  energy  terms,  and  by  adding  the  non- 
linearities  (up  to  second  order)  in  extension  and  torsion 
produced  by  bending,  ‘almost-exact’  beam  models  can  be 
developed.  For  larger  bending  deformations,  the  nonlin- 
earities  can  be  handled  by  using  the  ‘geometrically  exact’ 
beam  theories  [122],  or  within  a  multibody  type  formu¬ 
lation  [123]  which  allow  both  rigid  and  elastic  motion  of 
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structural  components. 

Since  the  mid  1980s,  there  was  an  emphasis  on  im¬ 
proving  cross-section  analyses  to  account  for  anisotropy 
produced  by  composite  materials.  Research  in  compos¬ 
ite  blade  modeling  has  been  reviewed  in  detail  by  several 
researchers,  most  recently  by  Friedmann  in  2004  [124]. 
The  basic  idea  has  been  to  perform  a  linear  or  non-linear 
two-dimensional  cross-sectional  analysis,  uncoupled  from 
the  ID  non-linear  beam  theories  mentioned  above,  with 
the  goal  of  obtaining  the  cross-sectional  elastic  constants 
and  warping.  It  has  been  shown  using  asymptotic  meth¬ 
ods  that  such  uncoupling  assumptions  hold  true  in  most 
cases  for  rotor  blades  [125].  The  cross-sectional  elastic 
constants  defining  the  stress-strain  relation  are  however 
fully  coupled.  For  low  frequency  dynamic  behavior,  a 
classical  analysis  with  a  single  shear  strain  measure  (as 
in  the  case  of  isotropic  beams)  can  be  reasonably  accu¬ 
rate,  as  long  as  the  coupled  elastic  constants  are  calcu¬ 
lated  properly  [126].  For  higher  accuracy,  as  required  for 
rotor  analysis,  the  strain  vector  can  include  two  addi¬ 
tional  strains  to  account  for  the  transverse  shears,  or  one 
additional  strain  to  account  for  the  Vlasov  or  restrained 
warping.  The  accuracy  of  the  model  depends  more  on 
the  accuracy  of  the  elastic  constants  than  the  number  of 
strain  measures.  There  are  several  methods  for  extract¬ 
ing  these  constants.  They  are:  engineering  methods  that 
work  well  for  certain  type  of  cross-sections  [127],  asymp¬ 
totic  methods  which  yield  closed  form  solutions  for  thin 
walled  beams  [128,  129],  and  finite  element  methods.  The 
latter  can  be  based  either  on  St.  Venant’s  principle  [130] 
or  on  asymptotic  methods  [131].  The  latter  reference  pro¬ 
vides  a  generalized  treatment,  where  the  cross-sectional 
equations  and  the  one  dimensional  blade  equations  are 
rigorously  deduced  from  three-dimensional  elasticity  the¬ 
ory.  These  analyses  or  static  measurements  are  used  to¬ 
day  to  obtain  the  section  properties  required  to  operate  a 
beam-theory  based  comprehensive  code.  The  dynamics 
of  the  blade  is  then  calculated,  subject  to  proper  model¬ 
ing  of  the  boundary  conditions.  The  latter  is  a  challeng¬ 
ing  task  because  of  nonlinearities  associated  with  the  root 
end  of  the  blades,  the  prime  examples  of  which  are  lag 
dampers  and  control  systems. 

Gandhi  and  Chopra  [132]  incorporated  an  analytical 
elastomeric  lag  damper  model  in  a  rotor  comprehensive 
analysis.  Finite  element  based  damper  models  have  also 
been  developed,  such  as  by  Smith  et  al  [133].  In  general, 
the  difficulties  with  modeling  lag  dampers  are  numerous. 
Non-linearities  associated  with  temperature,  amplitude, 
and  frequency  of  these  structures  create  unknown  param¬ 
eters  and  reduce  prediction  accuracy  of  the  blade  models. 
Complex  boundary  conditions  can  be  handled  in  a  gener¬ 
alized  manner  using  a  multibody  formulation  [134].  Orig¬ 
inally  developed  for  spacecraft  dynamics,  it  provides  a 
framework  to  treat  rigid  and  flexible  structural  elements 


undergoing  arbitrary  rotations  and  translations  with  re¬ 
spect  to  one  another.  In  addition,  unlike  conventional 
formulations,  multibody  formulations  enable  the  increase 
in  scope  of  analysis  without  having  to  re-derive  and  re¬ 
validate  the  equations  with  each  additional  feature.  One 
of  the  early  rotor  aeromechanical  stability  analysis  devel¬ 
oped  to  represent  arbitrary  rigid  and  flexible  dynamics 
was  GRASP  by  Hodges  et  al.  [135].  This  multibody  like 
analysis  was  designed  to  treat  modern  rotors  with  ar¬ 
bitrary  boundary  conditions  in  steady,  axial  flight  and 
ground  contact  conditions.  Multibody  formulations  en¬ 
able  large  deformation  problems  to  be  accurately  treated 
while  still  retaining  second-order  nonlinear  beam  theory 
(i.e.,  without  using  geometrically  exact  theories).  This  is 
done  by  simply  breaking  up  the  rotor  blade  into  multi¬ 
ple  bodies  each  undergoing  only  moderate  deformations 
in  its  local  frame.  The  total  deformation  is  obtained  by 
adding  the  local  deformations  to  the  end  contribution  of 
the  adjoining  parent  element  as  a  rigid  motion. 

DYMORE  [136]  and  MBDyn  [137]  are  examples  of 
modern  multibody  structural  codes.  These  codes  are  of¬ 
ten  geared  to  solve  the  full  non-linear  finite  element  equa¬ 
tions  using  time- marching  schemes.  This  is  due  to  the 
presence  of  algebraic  constraint  equations.  Time  march¬ 
ing  schemes  are  less  suitable  for  rotor  trim.  Artificial 
damping  is  often  required  to  suppress  the  natural  re¬ 
sponse  of  the  lowly  damped  modes,  particularly  the  lag 
mode.  The  damping  needs  to  be  subsequently  removed 
as  the  solution  progresses  to  arrive  at  the  correct  steady 
state  periodic  response.  In  addition,  the  use  of  algebraic 
equations  prevents  the  use  of  state  space  based  aeroelas- 
tic  stability  analysis.  The  accuracy  of  blade  loads  predic¬ 
tion,  however,  have  not  shown  any  significant  improve¬ 
ment,  so  far,  with  the  expansion  of  the  scope  of  structural 
modeling.  The  importance  of  detailed  boundary  condi¬ 
tion  modeling,  and  their  impact  on  the  current  deficien¬ 
cies  in  blade  loads  prediction,  remains  to  be  understood. 

The  development  of  coupled  rotor-fuselage  dynamic 
analysis  can  be  traced  back  to  the  early  work  of  Gersten- 
berger  and  Wood  [138]  in  1963.  The  work  was  based  on 
a  method  of  impedance  matching,  which  was  later  used 
extensively  during  the  1970s  and  80s  [139,  140,  141,  142]. 
Beginning  1984,  NASA  Langley  Research  Center  carried 
out  a  successful  Design  Analysis  Methods  for  Vibrations 
(DAMVIBS)  program,  with  the  active  participation  of 
four  major  rotorcraft  manufacturers  (Bell,  Boeing,  the 
former  McDonnell-Douglas,  and  Sikorsky).  As  part  of 
the  program  detailed  FEM  of  several  metal  and  compos¬ 
ite  helicopter  fuselages  were  constructed,  e.g.  the  AH- 
1G,  Model  360,  AH-64A,  OH-6A,  CH-47D,  S75,  D292, 
and  the  UH-60A.  An  overview  of  the  program  can  be 
found  in  Kvaternik  [143].  The  description  of  the  airframe 
FEM  models  as  developed  by  each  of  the  companies  can 
be  found  in  Refs.  [144,  145,  146,  147].  The  prevalent 
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methods  of  rotor-fuselage  dynamic  analysis  were  sum¬ 
marized  by  Kvaternik  et  al.  [148]. 

During  the  1990s,  iterative  or  fully  coupled  formula¬ 
tions  began  to  be  applied  to  the  rotor-fuselage  problem. 
An  early  work  was  that  of  Stephens  and  Peters  [149],  fol¬ 
lowed  by  Vellaichamy  and  Chopra  [150],  and  Chiu  and 
Friedmann  [151].  During  the  2000s,  detailed  3-D  fuse¬ 
lage  FEM  models  began  to  be  coupled  to  rotor  dynam¬ 
ics  using  these  formulations.  Cribbs  et  al.  [152]  studied 
active  vibration  control  techniques  in  the  fixed  frame  us¬ 
ing  a  MBB  BO- 105  type  four-bladed  hingeless  rotor  he¬ 
licopter.  The  DAMVIBS  AH-1G  fuselage  was  used  by 
Yeo  and  Chopra  [153]  to  carry  out  vibration  prediction 
studies  (see  Figs.  6  and  7).  A  similar  (to  DAMVIBS) 
SH-60  fuselage  was  later  used  by  Yang  et  al.  [154],  to 
study  the  vibration  impact  of  damaged  rotors.  Bauchau 
et  al.  [155]  applied  the  method  of  component  mode  syn¬ 
thesis  to  couple  fuselage  dynamics  to  a  multibody  rotor 
model. 

Except  for  the  two-bladed  teetering  AH-1G  heli¬ 
copter,  the  inclusion  of  fuselage  dynamics  had  a  minor 
effect  on  the  blade  loads.  Predicted  vibration  in  the  fuse¬ 
lage,  of  course,  relies  on  fuselage  dynamics.  As  mentioned 
earlier,  all  of  these  studies  focused  on  vibratory  forcing 
transmitted  via  the  shaft.  The  effect  of  rotor  wake  im¬ 
pingement  were  not  included. 

Aeroelastic  Phenomena  and  CFD/CSD  Predic¬ 
tions 

In  rotary  wing  aeroelasticity,  a  first  principles  solu¬ 
tion  requires  more  than  CFD  and  CSD.  It  requires,  in 
addition,  a  method  for  simultaneously  predicting  the  un¬ 
steady  control  pitch  angles  and  aircraft  operating  states. 
This  is  referred  to  as  the  Vehicle  Flight  Dynamics  (VFD) 
problem  in  this  paper.  The  requirement  arises  be¬ 
cause  unlike  fixed- wings  and  turbomachines,  the  dynamic 
loads,  flutter  (flap-lag),  and  air-resonance  characteristics 
of  the  rotor  blades  are  coupled  inherently  and  nonlin- 
early  to  vehicle  dynamics  via  the  control  pitch  angles 
prescribed  at  the  blade  root. 

In  steady  flight,  straight-level  or  turns,  VFD  deter¬ 
mines  the  operating  state  which  holds  the  entire  aircraft 
in  equilibrium.  The  operating  state  is  defined  by  the 
rotor/s  pitch  control  angles  (only  collective  in  the  case 
of  tail  rotors),  aircraft  attitude  angles  and  rates,  fuselage 
properties  (aerodynamic  interaction,  weight-distribution, 
tail  properties  and  setting,  compounding  and  special  con¬ 
figurations).  The  word  ‘coupled  trim’  is  used  by  dynam- 
icists  to  describe  the  VFD  associated  with  steady  flight. 
During  unsteady  flight,  VFD  is  a  generalization  of  cou¬ 
pled  trim,  and  determines  the  operating  state  which  fol¬ 
lows  an  intended  flight  dynamics  of  the  aircraft.  For  a 
general  maneuver,  these  intended  targets  can  be  a  pre¬ 


scribed  trajectory  in  space,  or  a  prescribed  mean  loading 
schedule  (e.g.  g-level  and  a  prescribed  set  of  hub  roll 
and  pitch  moment  variation  with  time).  The  optimal 
control  problem  of  determining  the  control  angles  based 
on  these  intended  targets  is  still  beyond  state-of-the-art, 
even  with  lower  order  free  wake  based  lifting-line  aero¬ 
dynamic  models. 

For  steady  flight,  the  VFD  problem  is  relatively  sim¬ 
ple.  The  trajectory  is  trivial,  and  the  mean  loading 
schedule  can  be  calculated  from  first  principles  using  a  six 
DOF  equilibrium  in  level  flight,  or  an  eight  DOF  equilib¬ 
rium  in  a  most  general  steady  maneuver  (a  coordinated 
helical  turn).  Calculating  VFD  (or  ‘coupled  trim’)  is  the 
first  step  for  predicting  performance,  loads,  vibration, 
stability,  and  acoustics.  An  innovative  method  for  si¬ 
multaneous  determination  of  VFD  during  CFD/CSD  for 
rotors  was  described  by  Tung  et  al.  in  1986  [156].  Re¬ 
ferred  to  in  this  paper  as  the  ‘delta  coupling’,  it  has  re¬ 
mained  the  most  efficient  approach  and  has  been  the  cor¬ 
nerstone  of  all  current  advances  in  rotorcraft  CFD/CSD. 
The  method  is  also  referred  to  as  ‘loose  coupling’,  how¬ 
ever,  it  bears  no  resemblance  with  CFD/CSD  loose  cou¬ 
pling  in  fixed  wings,  and  the  confusion  is  best  avoided. 

For  unsteady  maneuvers,  the  analysis  becomes  con¬ 
ceptually  simpler  when  the  rotor  control  angles  and  ve¬ 
hicle  dynamics  are  known.  The  VFD  is  then  uncoupled 
from  CFD/CSD  simulation.  Prescribing  VFD,  either 
from  measured  data  or  from  isolated  lower  order  pre¬ 
dictions,  reduces  the  CFD/CSD  problem  to  a  straight¬ 
forward,  partitioned,  numerical  integration  problem, 
similar  to  those  applied  in  fixed  wing  and  turbo  machin¬ 
ery  disciplines.  Except  that,  in  rotorcraft,  the  problem 
is  currently  solved  at  a  much  lower  fidelity  compared  to 
the  other  disciplines.  First,  the  structural  model  for  an 
isolated  rotor  blade  is  almost  always  a  beam,  compared 
to  3-D  structures  in  other  disciplines.  Second,  the  em¬ 
phasis  in  rotors  has  been  on  blade  loads,  not  flutter,  and 
therefore  errors  in  approximate  fluid-structure  coupling 
methods  and  low-order  (below  second  order)  temporal  in¬ 
tegration  have  been  of  lesser  consequence.  There  is  only 
a  single  example  of  CFD/CSD  coupling  performed  so  far 
(see  later)  for  a  realistic  unsteady  maneuver. 

The  application  of  delta  coupling  during  the  1990s 
yielded  limited  payoffs  compared  to  lifting-line  meth¬ 
ods  [157,  158,  159].  A  comprehensive  report  on  loads 
prediction  for  the  Research  Puma  helicopter  using  U.S. 
and  European  CFD/CSD  methods  of  the  time  is  docu¬ 
mented  in  Bousman  et  al.  [160]  in  1996.  Philippe  [161], 
in  1994,  presented  a  summary  of  other  European  efforts. 
The  key  drawback  at  the  time  was  the  lack  of  consis¬ 
tent  coupling  due  to  the  inability  of  the  CFD  solvers  in 
resolving  the  grid  and  compressibility  effects.  Renewed 
interest  in  the  2000s  was  spurred  by  the  emergence  of 
reliable  unsteady  Euler /RANS  solvers. 
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The  resurgence  of  rotary  wing  CFD/CSD/VFD,  in 
the  U.S.,  stems  from  the  systematic  work  of  Datta  et 
al.  [162,  163]  and  Potsdam  et  al.  [164].  The  advances  were 
enabled  in  substantial  part  by  the  reliable  and  repeat- 
able  flight  test  data  provided  by  the  NAS  A/ Army  UH- 
60A  Airloads  Flight  Test  Program.  Reference  [162]  dis¬ 
sected  the  coupled  aeroelastic  loads  problem  into  aerody¬ 
namic,  structural  dynamic,  and  trim  components;  identi¬ 
fied  the  key  mechanisms  of  rotor  vibration  at  high  speed, 
isolated  the  effects  of  wake  and  transonic  pitching  mo¬ 
ments,  and  clarified  the  key  contribution  of  Euler /RANS 
over  lifting-line  models.  Reference  [163]  brought  the 
components  back  together  again  and  demonstrated  im¬ 
proved  airloads  and  structural  loads  prediction  using 
CFD/Comprehensive  Analysis  coupling.  The  predictions 
were  compared  systematically  between  lifting-line  air¬ 
loads,  CFD  airloads,  and  measured  airloads.  The  CFD 
airloads  were  computed  using  near  blade  RANS  coupled 
to  free-wake  inflow  using  the  field  velocity  approach  by 
Sitaraman  and  Baeder  [165]. 

Potsdam  et  al.  [164]  focused  on  airloads.  However, 
the  work  extended  RANS  capabilities  beyond  the  high 
speed  flight  to  a  wake  dominated  transition  flight  and  a 
stall  dominated  highly  loaded  flight.  In  the  wake  dom¬ 
inated  flight,  the  work  demonstrated  that  CFD  is  capa¬ 
ble  of  predicting  the  dominant  wake  induced  vibratory 
airloads  at  transition  speeds  even  without  fully  captur¬ 
ing  the  individual  vortices  in  the  wake.  At  the  dynamic 
stall  flight,  the  work  demonstrated  that  a  one  equa¬ 
tion  turbulence  model  (Baldwin-Barth)  is  adequate  for 
predicting  the  basic  stall  phenomena  on  the  retreating 
side.  The  stall  flight  was  subsequently  studied  in  de¬ 
tail  by  Datta  and  Chopra  [166],  isolating  the  physics  of 
structural  dynamics,  wake,  stall,  and  trim.  The  physics 
of  the  stall  mechanisms  and  their  impact  on  predicted 
pitch-link  loads  using  CFD / Comprehensive  Analysis  cou¬ 
pling  were  clarified.  The  CFD  analysis  used  in  refer¬ 
ences  [162,  163,  166]  was  subsequently  extended  to  in¬ 
clude  overset  capabilities  for  wake  capture,  similar  to 
that  of  reference  [164],  by  Duraisamy  and  Baeder  [167]. 
This  methodology  was  then  used  to  carry  out  predictions 
of  vibratory  airloads  and  structural  loads  at  all  three  level 
flight  conditions  in  Refs.  [168,  169].  The  work  of  Pots¬ 
dam  et  al.  [164]  was  subsequently  extended  by  Lim  et 
al.  [170]  for  studying  vortex  loadings  in  more  details.  The 
focus  was  on  descent  flights  of  the  UH-60A  and  HART 
II  rotor  where,  both  the  vortex  induced  loadings  as  well 
as  the  individual  BVIs  are  present  in  the  airloads.  The 
BVI  airloads  predicted  by  CFD  were  improved  compared 
to  lifting-line  predictions.  A  significant  portion  of  the 
improvement  was  in  low  frequency  harmonics  (3/rev), 
brought  about  by  vortex  loadings  in  the  first  and  fourth 
quadrants,  rather  than  individual  BVI.  Increasing  the 
background  grid  resolution  from  10%  chord  to  5%  chord, 


a  difference  of  19. 4M  and  107. 6M  grid  points,  did  not 
show  any  improvement  in  the  individual  BVI  loadings. 

Concurrent  developments  in  Europe  over  the  last 
five  years  have  also  contributed  to  the  increased  un¬ 
derstanding  of  CFD/Comprehensive  Analysis  capabili¬ 
ties  for  loads  prediction  in  helicopter  rotors.  The  high 
speed  wind  tunnel  test  airloads  of  the  ONERA  7  A  and 
7AD  rotors  have  been  studied  and  predicted  by  sev¬ 
eral  researchers,  see  Servera  et  al.  [171],  Altmikus  et 
al  [172],  Pomin  and  Wagner  [173],  and  Pahlke  and  Van 
Der  Wall  [174].  A  discussion  on  the  key  accomplishe- 
ments  in  all  of  the  above  work  can  be  found  in  the  review 
article,  Ref.  [116]. 

The  developments  cited  above  provided  a  renewed 
impetus  for  integrating  CFD  reliably  within  the  rotor 
design  process.  At  the  time,  among  all  key  design  is¬ 
sues,  noise,  tied  directly  to  vulnerability,  was  considered 
to  be  of  immediate  urgency.  Innovative  blade  shapes 
and  advanced  blade  geometries  for  noise  reduction  are 
not  considered  as  serious  contenders  for  blade  design  be¬ 
cause  their  impact  on  performance,  loads,  and  vibration 
cannot  be  predicted  with  confidence.  As  CFD  begins  to 
predict  the  fundamental  mechanisms  for  structural  loads 
-  3-D  unsteady  transonic  pitching  moments,  wake  in¬ 
duced  airloads  and  BVI,  and  dynamic  stall  cycles  —  a 
considerable  attention  is  now  focused  on  increasing  CFD 
throughput  while  maintaining  the  same  level  of  accu¬ 
racy  as  benchmarked  by  the  above  cited  works.  The 
DARPA  Helicopter  Quieting  Program  was  designed  to 
carry  out  such  a  task.  A  key  objective  of  the  DARPA 
program  was  to  adapt  non-rotorcraft  CFD  with  demon¬ 
strated  scalability  to  rotor  airloads  prediction.  For  exam¬ 
ple,  the  TFLOW2000  code  was  adapted  for  rotary  wing 
CFD  (extended  and  renamed  SUmb)  and  coupled  to  a 
rotor  craft  comprehensive  analysis,  as  part  of  the  pro¬ 
gram.  An  additional  objective  of  the  program  was  to 
investigate  the  impact  of  high  fidelity  wake  and  turbu¬ 
lence  modeling,  e.g.  methods  for  solution  adaptive  back¬ 
ground  mesh  refinements,  different  near-blade  and  back¬ 
ground  solvers,  and  v2-f  and  KES  turbulence  models, 
which  these  analyses  incorporated.  The  focus  was  not  on 
CSD  or  coupling  methods.  To  this  end,  trimmed  aeroe¬ 
lastic  blade  deformations  from  Ref.  [164],  were  supplied 
to  the  participants  as  an  initial  step.  As  a  follow  on,  the 
participants  performed  the  same  CFD/CSD/VFD  cou¬ 
pling  as  in  Refs.  [163,  164]  using  the  same  comprehen¬ 
sive  analysis  tools.  The  validation  set  was  extended  to 
include  acoustic  signatures  in  addition  to  performance 
and  airloads.  The  DNW  model  UH-60A,  and  the  BO- 
105  model  HART-II  rotor  were  considered  as  valida¬ 
tion  cases  for  acoustics.  The  contributions  of  the  three 
program  participant  teams  are  partially  documented  in 
Refs.  [176,  177,  178,  179]. 

All  of  the  above  coupled  studies,  with  the  exception 
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of  reference  [173],  used  variations  of  the  delta  coupling 
method  in  the  time  domain.  Reference  [173]  does  not 
calculate  control  angles  as  part  of  the  CFD  coupled  so¬ 
lution  but  determines  them  a  priori  using  an  uncoupled 
lifting-line  comprehensive  analysis.  The  delta  method  is 
unique  to  rotorcraft  and  is  not  found  in  fixed-wing  or 
turbomachinary  aeroelasticity.  It  accounts  for  the  fol¬ 
lowing  key  requirements  of  rotor  CSD:  (1)  it  operates, 
by  design,  at  or  near  resonance  in  flap,  (2)  it  requires 
simultaneous  solution  of  VFD  for  the  pitch  control  an¬ 
gles.  The  description  of  the  original  method  can  be 
found  in  Ref.  [156],  its  current  adaptations  in  any  of 
the  above  references,  see  for  example  Refs.  [164,  163]. 
The  method  is  implied  in  rotorcraft  literature  by  the  fre¬ 
quent  use  of  the  term  ‘CFD/Comprehensive  Analysis’. 
The  term  ‘Comprehensive  Analysis’,  by  itself,  histori¬ 
cally  refers  to  a  complete  analysis,  which  contains  within 
itself  all  of  the  airloads,  structural  loads,  and  VFD  mod¬ 
els.  Thus,  ‘CFD/Comprehensive  Analysis’  simply  refers 
to  a  CFD/CSD/VFD  analysis  which  uses  the  CSD  and 
VFD  of  existing  comprehensive  analysis  while  replacing 
the  built-in  lifting-line/surface  aerodynamics  with  CFD 
via  the  delta  coupling  method. 

In  the  special  case  where  the  control  angle  varia¬ 
tion  is  assumed  to  be  known  a  priori,  either  from  mea¬ 
sured  test  data  or  from  separate  uncoupled  estimates 
from  lower  order  models,  the  problem  reduces  to  a  pure 
CFD/CSD  analysis.  A  straight-forward,  partitioned,  nu¬ 
merical  integration  procedure  can  then  be  applied,  sim¬ 
ilar  to  those  used  in  fixed-wing  or  turbomachinery.  Ex¬ 
cept  that,  for  an  isolated  rotor  analysis,  the  problem  is 
rendered  much  simpler  as  the  structural  model  is  almost 
always  a  beam,  and  a  low  order  fluid-structure  interface 
is  much  simpler  to  devise.  These  methods  are  known  was 
‘tight  coupling’  in  rotorcraft.  Again,  note  the  confusing 
nomenclature  as  compared  to  fixed  wing  terminology.  In 
fixed  wing  tight  coupling  implies  strict  time  accuracy  via 
sub- iterations.  One  of  the  first  implementations  of  tight 
coupling  in  rotorcraft  using  RANS  CFD  was  by  Bauchau 
and  Ahmad  [180].  Recent  studies  are  those  of  Pomin  and 
Wagner  [173]  and  Bhagwat  et  al  [181].  The  latter  work 
is  of  significance  as  it  is  the  first  application  of  CFD  to 
a  real  unsteady  maneuver.  The  method  is  applied  here 
to  a  pull-up  maneuver,  the  second  most  severe  (loads) 
of  the  UH-60A  maneuvers.  Note  that,  the  advances  in 
this  study  were  enabled  by  the  NASA/Army  UH-60A 
Airloads  Flight  Test  Program  Data.  The  rotor  control 
angles  and  vehicle  dynamics  were  prescribed  from  mea¬ 
sured  flight  data.  Several  key  conclusions  can  be  drawn 
from  this  work.  First,  it  demonstrates  RANS  capability 
for  predicting  advancing  blade  transonic  stall.  Second, 
it  appears  to  negate  the  role  of  fuselage  upwash  on  ro¬ 
tor  stall  and  pitch-link  loads  during  a  pull-up  maneuver. 
Third,  it  demonstrates  that  an  isolated  rotor  calculation 


can  predict  the  oscillatory  (peak-to-peak)  control  loads 
correctly  in  the  severest  of  the  UH-60A  maneuvers. 

A  LIST  OF  TEMPORAL  COUPLING 
PROCEDURES 

In  view  of  the  conflicting  nomenclatures  often  used 
across  disciplines  to  describe  the  various  methods  of  tem¬ 
poral  coupling,  a  list  is  summarized  below. 

1.  Loose  coupling :  Partitioned  formulation  based  tem¬ 
poral  advancement,  but  without  sub-iterations. 
Predictor-corrector  schemes  often  used  to  seek  near- 
second  order  accuracy.  Standard  nomenclature  in 
fixed  wing.  Described  as  tight  coupling  in  rotary 
wing. 

2.  Tight  coupling :  Partitioned  formulation  based  time 
advancement,  but  with  sub- iterations.  Temporal  ac¬ 
curacy  level  of  worse  solver.  Standard  nomenclature 
in  fixed  wing.  Also  called  close  coupling  or  strong 
coupling. 

3.  Periodic  coupling :  For  periodic  dynamics.  First  ap¬ 
plied  by  Gerolymos  for  turbomachinery.  Iterates  on 
harmonics  of  loads  and  deformations.  Strict  time 
accuracy  can  be  enforced  on  any  chosen  set  of  har¬ 
monics. 

4.  Delta  coupling :  Periodic  coupling  with  the  use  of 
aerodynamic  damping  from  lower  order  models. 
Only  successful  approach  for  CFD/CSD/VFD  in  ro¬ 
torcraft.  Conflictingly  nomenclature  as  loose  cou¬ 
pling  in  rotary  wing.  Note,  not  loose  in  the  fixed 
wing  sense  of  the  term.  Strict  time  accuracy  can  be 
enforced  on  any  chosen  set  of  harmonics. 

5.  Full  coupling :  Formulation  coupled  via  boundary 
conditions  (deflection  and  velocity  compatibility, 
and  stress  equilibrium)  and  integrated  together. 
Also  called  direct  coupling ,  intimate  coupling ,  or 
monolithic  coupling. 

CONCLUDING  OBSERVATIONS 

The  objective  of  the  present  paper  was  to  moti¬ 
vate  research  towards  the  next  generation  HPC  based 
high-fidelity  rotorcraft  comprehensive  analysis.  To  this 
end,  the  frontiers  of  fixed  wing,  turbomachinery,  and  ro¬ 
tary  wing  aeroelasticity  were  reviewed,  along  with  a  sur¬ 
vey  of  the  fluid-structure  coupling  methodologies  used 
in  each.  The  underlying  theme  of  the  paper  was  rotor¬ 
craft.  Therefore,  the  technology  areas  in  each  discipline 
were  tied  to  the  future  rotorcraft  requirements,  where 
deemed  relevant  by  the  authors.  In  addition,  compar¬ 
isons  were  drawn  between  the  computational  procedures 
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used  in  each  of  the  disciplines,  and  those  required  by  the 
unique  rot  or  craft  requirements. 

CFD  is  only  one  part  of  a  next  generation  com¬ 
prehensive  analysis.  The  ultimate  objective  of  reliably 
discriminating  between  innovative  blade  shapes,  rotor 
types,  and  rotor-fuselage  configurations,  requires  CFD, 
CSD,  high-fidelity  coupling,  and  a  solution  procedure 
that  ties  CFD/CSD  to  VFD  satisfying  the  unique  oper¬ 
ating  principles  of  rotorcraft.  The  focus  so  far  has  largely 
been  on  CFD.  This  is  because  CFD,  up  to  recently,  was 
unable  to  meet  the  fundamental  requirements  for  vibra¬ 
tory  structural  loads  in  rotorcraft:  predicting  transonic 
pitching  moments  on  blades,  wake  induced  airloads  on 
both  blades  and  fuselage,  dynamic  stall  cycles  on  blades, 
and  buffet  loads  on  the  fuselage  driven  by  massive  sep¬ 
aration.  HPC  based  CFD  is  now  beginning  to  address 
these  requirements.  The  methods  of  high-fidelity  CSD, 
fluid-structure  coupling,  and  solution  procedures  must 
now  be  re-evaluated  and  re-shaped  for  taking  rotorcraft 
comprehensive  analysis  to  the  next  level. 
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